{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b8bba8eb",
   "metadata": {},
   "source": [
    "# Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cddf1967",
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import basename, exists\n",
    "\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print(\"Downloaded \" + local)\n",
    "\n",
    "\n",
    "download(\"https://github.com/AllenDowney/ThinkStats2/raw/master/code/thinkstats2.py\")\n",
    "download(\"https://github.com/AllenDowney/ThinkStats2/raw/master/code/thinkplot.py\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bb3a6450",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import thinkstats2\n",
    "import thinkplot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dacf078",
   "metadata": {},
   "source": [
    "## The estimation game\n",
    "\n",
    "Let's play a game. I think of a distribution, and you have to guess what\n",
    "it is. I'll give you two hints: it's a normal distribution, and here's a\n",
    "random sample drawn from it:\n",
    "\n",
    "`[-0.441, 1.774, -0.101, -1.138, 2.975, -2.138]`\n",
    "\n",
    "What do you think is the mean parameter, $\\mu$, of this distribution?\n",
    "\n",
    "One choice is to use the sample mean, $\\bar{x}$, as an estimate of $\\mu$.\n",
    "In this example, $\\bar{x}$ is 0.155, so it would be reasonable to guess\n",
    "$\\mu$ = 0.155. This process is called **estimation**, and the statistic\n",
    "we used (the sample mean) is called an **estimator**.\n",
    "\n",
    "Using the sample mean to estimate $\\mu$ is so obvious that it is hard to\n",
    "imagine a reasonable alternative. But suppose we change the game by\n",
    "introducing outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f2a7144",
   "metadata": {},
   "source": [
    "*I'm thinking of a distribution.* It's a normal distribution, and here's\n",
    "a sample that was collected by an unreliable surveyor who occasionally\n",
    "puts the decimal point in the wrong place.\n",
    "\n",
    "`[-0.441, 1.774, -0.101, -1.138, 2.975, -213.8]`\n",
    "\n",
    "Now what's your estimate of $\\mu$? If you use the sample mean, your\n",
    "guess is -35.12. Is that the best choice? What are the alternatives?\n",
    "\n",
    "One option is to identify and discard outliers, then compute the sample\n",
    "mean of the rest. Another option is to use the median as an estimator.\n",
    "\n",
    "Which estimator is best depends on the circumstances (for example,\n",
    "whether there are outliers) and on what the goal is. Are you trying to\n",
    "minimize errors, or maximize your chance of getting the right answer?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f619e94",
   "metadata": {},
   "source": [
    "If there are no outliers, the sample mean minimizes the **mean squared\n",
    "error** (MSE). That is, if we play the game many times, and each time\n",
    "compute the error $\\bar{x} - \\mu$, the sample mean minimizes\n",
    "$$MSE = \\frac{1}{m} \\sum (\\bar{x} - \\mu)^2$$ Where $m$ is the number of\n",
    "times you play the estimation game, not to be confused with $n$, which\n",
    "is the size of the sample used to compute $\\bar{x}$.\n",
    "\n",
    "Here is a function that simulates the estimation game and computes the\n",
    "root mean squared error (RMSE), which is the square root of MSE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c6ae3ee6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def RMSE(estimates, actual):\n",
    "    \"\"\"Computes the root mean squared error of a sequence of estimates.\n",
    "\n",
    "    estimate: sequence of numbers\n",
    "    actual: actual value\n",
    "\n",
    "    returns: float RMSE\n",
    "    \"\"\"\n",
    "    e2 = [(estimate-actual)**2 for estimate in estimates]\n",
    "    mse = np.mean(e2)\n",
    "    return np.sqrt(mse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7c03dd73",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "def Estimate1(n=7, m=1000):\n",
    "    mu = 0\n",
    "    sigma = 1\n",
    "\n",
    "    means = []\n",
    "    medians = []\n",
    "    for _ in range(m):\n",
    "        xs = [random.gauss(mu, sigma) for i in range(n)]\n",
    "        xbar = np.mean(xs)\n",
    "        median = np.median(xs)\n",
    "        means.append(xbar)\n",
    "        medians.append(median)\n",
    "\n",
    "    print('rmse xbar', RMSE(means, mu))\n",
    "    print('rmse median', RMSE(medians, mu))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6484add8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse xbar 0.39530236044347977\n",
      "rmse median 0.475238334720722\n"
     ]
    }
   ],
   "source": [
    "Estimate1()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9700bb61",
   "metadata": {},
   "source": [
    "Again, `n` is the size of the sample, and `m` is the number of times we\n",
    "play the game. `means` is the list of estimates based on $\\bar{x}$.\n",
    "`medians` is the list of medians.\n",
    "\n",
    "Here's the function that computes RMSE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f4a0e261",
   "metadata": {},
   "outputs": [],
   "source": [
    "def RMSE(estimates, actual):\n",
    "    e2 = [(estimate-actual)**2 for estimate in estimates]\n",
    "    mse = np.mean(e2)\n",
    "    return np.sqrt(mse)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf65739b",
   "metadata": {},
   "source": [
    "`estimates` is a list of estimates; `actual` is the actual value being\n",
    "estimated. In practice, of course, we don't know `actual`; if we did, we\n",
    "wouldn't have to estimate it. The purpose of this experiment is to\n",
    "compare the performance of the two estimators.\n",
    "\n",
    "When I ran this code, the RMSE of the sample mean was 0.41, which means\n",
    "that if we use $\\bar{x}$ to estimate the mean of this distribution, based\n",
    "on a sample with $n=7$, we should expect to be off by 0.41 on average.\n",
    "Using the median to estimate the mean yields RMSE 0.53, which confirms\n",
    "that $\\bar{x}$ yields lower RMSE, at least for this example.\n",
    "\n",
    "Minimizing MSE is a nice property, but it's not always the best\n",
    "strategy. For example, suppose we are estimating the distribution of\n",
    "wind speeds at a building site. If the estimate is too high, we might\n",
    "overbuild the structure, increasing its cost. But if it's too low, the\n",
    "building might collapse. Because cost as a function of error is not\n",
    "symmetric, minimizing MSE is not the best strategy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d994ac2",
   "metadata": {},
   "source": [
    "As another example, suppose I roll three six-sided dice and ask you to\n",
    "predict the total. If you get it exactly right, you get a prize;\n",
    "otherwise you get nothing. In this case the value that minimizes MSE is\n",
    "10.5, but that would be a bad guess, because the total of three dice is\n",
    "never 10.5. For this game, you want an estimator that has the highest\n",
    "chance of being right, which is a **maximum likelihood estimator**\n",
    "(MLE). If you pick 10 or 11, your chance of winning is 1 in 8, and\n",
    "that's the best you can do."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "981535aa",
   "metadata": {},
   "source": [
    "## Guess the variance\n",
    "\n",
    "*I'm thinking of a distribution.* It's a normal distribution, and here's\n",
    "a (familiar) sample:\n",
    "\n",
    "`[-0.441, 1.774, -0.101, -1.138, 2.975, -2.138]`\n",
    "\n",
    "What do you think is the variance, $\\sigma^2$, of my distribution?\n",
    "Again, the obvious choice is to use the sample variance, $S^2$, as an\n",
    "estimator. $$S^2 = \\frac{1}{n} \\sum (x_i - \\bar{x})^2$$ For large samples,\n",
    "$S^2$ is an adequate estimator, but for small samples it tends to be too\n",
    "low. Because of this unfortunate property, it is called a **biased**\n",
    "estimator. An estimator is **unbiased** if the expected total (or mean)\n",
    "error, after many iterations of the estimation game, is 0.\n",
    "\n",
    "Fortunately, there is another simple statistic that is an unbiased\n",
    "estimator of $\\sigma^2$:\n",
    "$$S_{n-1}^2 = \\frac{1}{n-1} \\sum (x_i - \\bar{x})^2$$ For an explanation of\n",
    "why $S^2$ is biased, and a proof that $S_{n-1}^2$ is unbiased, see\n",
    "<http://wikipedia.org/wiki/Bias_of_an_estimator>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17f227d8",
   "metadata": {},
   "source": [
    "The biggest problem with this estimator is that its name and symbol are\n",
    "used inconsistently. The name \"sample variance\" can refer to either\n",
    "$S^2$ or $S_{n-1}^2$, and the symbol $S^2$ is used for either or both.\n",
    "\n",
    "Here is a function that simulates the estimation game and tests the\n",
    "performance of $S^2$ and $S_{n-1}^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "48642d63",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MeanError(estimates, actual):\n",
    "    \"\"\"Computes the mean error of a sequence of estimates.\n",
    "\n",
    "    estimate: sequence of numbers\n",
    "    actual: actual value\n",
    "\n",
    "    returns: float mean error\n",
    "    \"\"\"\n",
    "    errors = [estimate-actual for estimate in estimates]\n",
    "    return np.mean(errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "09615c49",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Estimate2(n=7, m=1000):\n",
    "    mu = 0\n",
    "    sigma = 1\n",
    "\n",
    "    estimates1 = []\n",
    "    estimates2 = []\n",
    "    for _ in range(m):\n",
    "        xs = [random.gauss(mu, sigma) for i in range(n)]\n",
    "        biased = np.var(xs)\n",
    "        unbiased = np.var(xs, ddof=1)\n",
    "        estimates1.append(biased)\n",
    "        estimates2.append(unbiased)\n",
    "\n",
    "    print('mean error biased', MeanError(estimates1, sigma**2))\n",
    "    print('mean error unbiased', MeanError(estimates2, sigma**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c75bc848",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean error biased -0.13695189314208897\n",
      "mean error unbiased 0.006889458000896184\n"
     ]
    }
   ],
   "source": [
    "Estimate2()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27caa4cd",
   "metadata": {},
   "source": [
    "Again, `n` is the sample size and `m` is the number of times we play the\n",
    "game. `np.var` computes $S^2$ by default and $S_{n-1}^2$ if you provide\n",
    "the argument `ddof=1`, which stands for \"delta degrees of freedom.\" I\n",
    "won't explain that term, but you can read about it at\n",
    "<http://en.wikipedia.org/wiki/Degrees_of_freedom_(statistics)>.\n",
    "\n",
    "`MeanError` computes the mean difference between the estimates and the\n",
    "actual value:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b56aab89",
   "metadata": {},
   "source": [
    "When I ran this code, the mean error for $S^2$ was -0.13. As expected,\n",
    "this biased estimator tends to be too low. For $S_{n-1}^2$, the mean\n",
    "error was 0.014, about 10 times smaller. As `m` increases, we expect the\n",
    "mean error for $S_{n-1}^2$ to approach 0.\n",
    "\n",
    "Properties like MSE and bias are long-term expectations based on many\n",
    "iterations of the estimation game. By running simulations like the ones\n",
    "in this chapter, we can compare estimators and check whether they have\n",
    "desired properties.\n",
    "\n",
    "But when you apply an estimator to real data, you just get one estimate.\n",
    "It would not be meaningful to say that the estimate is unbiased; being\n",
    "unbiased is a property of the estimator, not the estimate.\n",
    "\n",
    "After you choose an estimator with appropriate properties, and use it to\n",
    "generate an estimate, the next step is to characterize the uncertainty\n",
    "of the estimate, which is the topic of the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ae4b7f1",
   "metadata": {},
   "source": [
    "## Sampling distributions\n",
    "\n",
    "Suppose you are a scientist studying gorillas in a wildlife preserve.\n",
    "You want to know the average weight of the adult female gorillas in the\n",
    "preserve. To weigh them, you have to tranquilize them, which is\n",
    "dangerous, expensive, and possibly harmful to the gorillas. But if it is\n",
    "important to obtain this information, it might be acceptable to weigh a\n",
    "sample of 9 gorillas. Let's assume that the population of the preserve\n",
    "is well known, so we can choose a representative sample of adult\n",
    "females. We could use the sample mean, $\\bar{x}$, to estimate the unknown\n",
    "population mean, $\\mu$.\n",
    "\n",
    "Having weighed 9 female gorillas, you might find $\\bar{x}=90$ kg and\n",
    "sample standard deviation, $S=7.5$ kg. The sample mean is an unbiased\n",
    "estimator of $\\mu$, and in the long run it minimizes MSE. So if you\n",
    "report a single estimate that summarizes the results, you would report\n",
    "90 kg.\n",
    "\n",
    "But how confident should you be in this estimate? If you only weigh\n",
    "$n=9$ gorillas out of a much larger population, you might be unlucky and\n",
    "choose the 9 heaviest gorillas (or the 9 lightest ones) just by chance.\n",
    "Variation in the estimate caused by random selection is called\n",
    "**sampling error**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24339e36",
   "metadata": {},
   "source": [
    "To quantify sampling error, we can simulate the sampling process with\n",
    "hypothetical values of $\\mu$ and $\\sigma$, and see how much $\\bar{x}$\n",
    "varies.\n",
    "\n",
    "Since we don't know the actual values of $\\mu$ and $\\sigma$ in the\n",
    "population, we'll use the estimates $\\bar{x}$ and $S$. So the question we\n",
    "answer is: \"If the actual values of $\\mu$ and $\\sigma$ were 90 kg and\n",
    "7.5 kg, and we ran the same experiment many times, how much would the\n",
    "estimated mean, $\\bar{x}$, vary?\"\n",
    "\n",
    "The following function answers that question:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "8f35a3ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def SimulateSample(mu=90, sigma=7.5, n=9, iters=1000):\n",
    "    xbars = []\n",
    "    for j in range(iters):\n",
    "        xs = np.random.normal(mu, sigma, n)\n",
    "        xbar = np.mean(xs)\n",
    "        xbars.append(xbar)\n",
    "    return xbars"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c43ec017",
   "metadata": {},
   "source": [
    "`mu` and `sigma` are the *hypothetical* values of the parameters. `n` is\n",
    "the sample size, the number of gorillas we measured. `m` is the number\n",
    "of times we run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2a26d949",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/K0lEQVR4nO3deZyP9f7/8ednhlkwi3UwhqEkWohwRttJU6NFOic1ISSk0slSJ1RMHdVUIi2KlFBZyjmpX6STiU6hFNE5LcoWxYx9Zsww6/v3R19T0+f6zMJnruuzPO632/zxeb2v6/q8PpfPXPN0rS5jjBEAAECACHG6AQAAAG8i3AAAgIBCuAEAAAGFcAMAAAIK4QYAAAQUwg0AAAgohBsAABBQajndgN1KS0u1Z88eRUVFyeVyOd0OAACoAmOMcnNz1bx5c4WEVLxvJujCzZ49e5SQkOB0GwAA4CTs3r1bLVq0qHCaoAs3UVFRkn5dOdHR0Q53AwAAqiInJ0cJCQllf8crEnTh5sShqOjoaMINAAB+piqnlHBCMQAACCiEGwAAEFAINwAAIKAQbgAAQEAh3AAAgIBCuAEAAAGFcAMAAAIK4QYAAAQUwg0AAAgohBsAABBQHA03//nPf9S7d281b95cLpdLS5curXSe1atXq3PnzgoPD9fpp5+uuXPn1nifAADAfzgabvLy8tSxY0fNmDGjStPv2LFDV199tS699FJt2rRJo0eP1rBhw/TBBx/UcKcAAKAieccKtDvzcNlP5oEcx3px9MGZV155pa688soqTz9z5ky1bt1aU6dOlSS1b99en376qZ5++mmlpKTUVJsAAPiM4uISZR3KdauXlhp98d+dKiwuVojFwyWLikr0r5VfqUmDKIWH1/Z6X7v3Hir3umWzBnp6/I1ef5+q8Kungq9bt07JycnlaikpKRo9erTHeQoKClRQUFD2OifHuSQJAAheOUeP6ZMNW3U0v0Chob8dONm195DWbNyqls0aVLqMXX8IECdjn0UwCjR+FW4yMzMVFxdXrhYXF6ecnBwdO3ZMkZGRbvOkp6fr4YcftqtFAAAk/RpmjuYXaOrcldr5y4FKp/dGcPEl0fUiHHtvvwo3J2PChAkaO3Zs2eucnBwlJCQ42BEAwJ8VF5foYHaeJKmwqEQbvvlJpaWmbHzrrn36/OsdTrXnE1ySrrzobMfe36/CTdOmTZWVlVWulpWVpejoaMu9NpIUHh6u8PBwO9oDAASovGMF+tuji5Sde8zpVqrsT+e2tqzvP3xUp7dsog6nN6uR93W5XDq9ZWPFNYyukeVXhV+Fm6SkJC1fvrxc7cMPP1RSUpJDHQEAAlVxcYmenp+hzzZvr5Hldz07sdzrPfuO6LKk9moYU7fSeevWCVf7Nk0VVjvUbczlcsllcUJxMHE03Bw9elRbt24te71jxw5t2rRJDRo0UMuWLTVhwgT98ssvmj9/viTp9ttv1/PPP6/77rtPt956qz766CO9+eabWrZsmVMfAQAQoFLvme31ZU4bd4NaNmsQ9OGjpjkabr788ktdeumlZa9PnBszePBgzZ07V3v37tWuXbvKxlu3bq1ly5ZpzJgxeuaZZ9SiRQu9/PLLXAYOAPCaj7/4Qc++/tFJzdulQytJkpHR7r2HddeAP+uMxDjVrhVKoLGRyxhjKp8scOTk5CgmJkbZ2dmKjnbueCAAwLes27RdT7367ypNO25YL511ejPV+r9LusNq1yK81LDq/P32q3NuAADwlpKSUj35ygc6lJOv3XsPqai4pMLp7+p/qbqek6h6dbhIxdcRbgAAQeN4QZHe+/i/envlVzpeUFTl+eY/PkR1Iwk1/oJwAwAIWD/+lKXJLy5X6xYN9b8f91R7/jGDk3Vh59NroDPUJMINAMDvGWOUf7yw7PWOnw8o7fn/V/a6usHm4vPb6u6be3IejZ8i3AAA/Nqyj/+rOf9ac8rLGfKXHqpdK1TndWipJg2ivNAZnEK4AQD4pa+3/KyHX3jvlJZx1cVnq13rpurcvqXqRIZ5qTM4jXADAPBpxcUl+nbbXh08kqf3P/mfjuTm6+CRvGovp+vZiTq7bXPVjQzXRV1OV61a7nf3RWAg3AAAfEppaam+256pnb8c9MrhpuSk9ro99WLOnwkihBsAgM948Jl39N32vSc9/7hhvdSmRSNJUnhYLUXVjfBWa/AjhBsAgGPyjxXqm217tOGbn/Th2u9Oejl/7tZOd/X/M3tnIIlwAwBwyN+f+qe2795f7fnq1QlXdN0IXdKtnZI6tVF8k1jvNwe/RrgBANiqpKRUN459qVrzXHFBB91yXZLCw2rXUFcIJIQbAIBtVnzyjWYv+aRK004bd4NaNW9Ywx0hEBFuAAA1rqioRDfdO7vCaUJDQ1RSUqoxg5J1QefTOH8GJ41wAwCoUVW52d6cRwYrJirSpo4Q6Ag3AIAa8dnm7Zoy59+VTrdgylDOpYFXEW4AAF738Rc/6NnXP6pwGg4/oaYQbgAAXpN/rFDvrNqsJR9sqHC6JdNHEGpQYwg3AACvWJqxSa+9+1mF0wzqk6Q+PTva1BGCFeEGAHDKPtu8vcJg0+vCszT8hots7AjBjHADADgpxhhNm7dSa7/aVuF0V150tob1vdCmrgDCDQDgJPy055DGPvFmpdONveVyXXDeaTZ0BPyGcAMAqJaRkxco80BOhdP0u7qb+l7R2aaOgPIINwCAKjHGqO/oWZVOR7CB0wg3AIAqqSzYjBvWS+edmaDatUNt6giwRrgBAFTq+lEzPY6lXHCWbruRK6HgOwg3AIAKVRRs3nhyqCLCeXQCfAvhBgDgUUXBhrsMw1eFON0AAMA3rang/jULnxpGsIHPItwAANwcLyjStLkfWo5NG3eDwmqz4x++i28nAMDNgPtesawvfGoYwQY+jz03AIByNn67y7J+x02XEGzgFwg3AIAyX323W4/OWm45dtmfzrS5G+DkEG4AAJKkwzn5emTmMssxroyCPyHcAAAkSZOefcey3vvP5xJs4FcINwAASdKe/dlutbatmuiWv/RwoBvg5BFuAABa8N56y3r6mL/Y3Alw6jjtHQCCnKe7ED80sjeHo+CX2HMDAEHq6y0/V/h4hXPOiLexG8B72HMDAEHonx9u9HgoSvp1rw3gr9hzAwBBpqiopMJgc1f/S9lrA7/GnhsACCIHjxzVbWmvexxfPHW4atUKtbEjwPsINwAQRDwFm0F9ktSnZ0ebuwFqBoelACBIrN20zeMYwQaBhHADAEGgtLRUU1/90HLszWm32dwNULMINwAQBO5+bLFlffHU4QoN5U8BAgvfaAAIAnstHq1wV/9LOXkYAYlwAwAB7tMNWy3rl3ZvZ3MngD0INwAQwIwxenr+Srf6Xf0vdaAbwB6EGwAIUIVFxeo7epblWOcOLW3uBrAP4QYAApAxRv3ufdlybMzgZMVERdrcEWAfwg0ABKCn52d4HLuw8+k2dgLYj3ADAAHmcE6+1my0Pol4/uNDbO4GsB+PXwCAAHK8oEjDJs63HFsyfYRcLpfNHQH2Y88NAASQAfe9Yll/YMRVBBsEDcINAAQIT/ezkbg6CsGFcAMAASDrYI7l/WykXw9HAcHE8XAzY8YMJSYmKiIiQt27d9f69esrnH769Olq166dIiMjlZCQoDFjxuj48eM2dQsAvifrYI7u/McCy7HFU4dzOApBx9Fws3jxYo0dO1ZpaWnauHGjOnbsqJSUFO3bt89y+gULFmj8+PFKS0vTd999p1deeUWLFy/W/fffb3PnAOA7PAWbQX2SeHYUgpKj4WbatGkaPny4hgwZog4dOmjmzJmqU6eO5syZYzn92rVrdcEFF6h///5KTEzUFVdcoX79+lW4t6egoEA5OTnlfgAgUOQcPWZZbxhbV316drS5G8A3OBZuCgsLtWHDBiUnJ//WTEiIkpOTtW7dOst5evTooQ0bNpSFme3bt2v58uW66qqrPL5Penq6YmJiyn4SEhK8+0EAwEFDHphnWX/p4YE2dwL4Dsfuc3PgwAGVlJQoLi6uXD0uLk7ff/+95Tz9+/fXgQMHdOGFF8oYo+LiYt1+++0VHpaaMGGCxo4dW/Y6JyeHgAMgIDy/YJVlffHU4TZ3AvgWx08oro7Vq1frscce0wsvvKCNGzfqX//6l5YtW6bJkyd7nCc8PFzR0dHlfgDAnx3KztMNo2dp1edb3Mbim8Ryng2CnmN7bho1aqTQ0FBlZWWVq2dlZalp06aW80ycOFEDBw7UsGHDJEnnnHOO8vLydNttt+mBBx5QSIhfZTUAqLZ9h3J1x8NveByfPuFGG7sBfJNjaSAsLExdunRRRsZvD3crLS1VRkaGkpKSLOfJz893CzChob/+D8UYU3PNAoAPMMZUGGxmpg3gP3mAHH621NixYzV48GCdf/756tatm6ZPn668vDwNGfLrg90GDRqk+Ph4paenS5J69+6tadOm6bzzzlP37t21detWTZw4Ub179y4LOQAQiEpLS3XDmJc8jl93WSc1bhBlY0eA73I03KSmpmr//v2aNGmSMjMz1alTJ61YsaLsJONdu3aV+1/Igw8+KJfLpQcffFC//PKLGjdurN69e+vRRx916iMAgC0qCjbzHx+iupHhNnYD+DaXCbLjOTk5OYqJiVF2djYnFwPwC0fzCzR4wquWY/985nabuwGcUZ2/3xycBQAf5ynYvJg2wOZOAP/g6GEpAMDJWTx1OJd8Ax6w5wYAfNje/dlutfCw2gQboAKEGwDwYXc9stCt9spkHq0AVIRwAwA+auHyLyzrkRFhNncC+BfCDQD4oB0/H9CSDzY43Qbglwg3AOCD7p2yxLL+xpNDbe4E8D+EGwDwMa+/+5ll/dkHblJEeG2buwH8D+EGAHzI7szDejtjk1u9x3mnKb5JrO39AP6IcAMAPmR0+mLL+j23XG5zJ4D/ItwAgI8oLCq2rD8w4iqbOwH8G+EGAHxEv3tfdqu1btFInTu0dKAbwH8RbgDAB8x+6xPL+pP3/NXmTgD/R7gBAIcZY7Ti028sx0JC2EwD1cVvDQA4bNuu/Zb1JdNH2NwJEBgINwDgsO+2Z7rVpo27US6Xy4FuAP9HuAEABxljNHfpWrd6q+YNHOgGCAyEGwBwUN/Rs9xq3IUYODWEGwBwyI6fD1jWB1zTzeZOgMBCuAEAh8xYuNqyftXF59jbCBBgCDcA4BCrPTeLpw53oBMgsBBuAMABVo9aaN44RrVqhTrQDRBYCDcA4IBvtu51q40aeJkDnQCBh3ADADbLzj2mR2Yuc6uf3qqJA90AgYdwAwA2u/XBeU63AAQ0wg0A2Oi5N1ZZ1mOj6tjcCRC4CDcAYJPjBUVavX6L5diLaf1t7gYIXLWcbgAAgsWA+16xrL857TaFhvJ/TcBb+G0CABscys6zrP/91isINoCX8RsFADYYPuk1y/qfOraxuRMg8BFuAKCG5R0rsKz/85nbbe4ECA6EGwCoYYPGv+p0C0BQIdwAQA26ftRMy/q89CE2dwIED8INANSQjM++s6zXqxOuenXCbe4GCB6EGwCoAUdy8/XCwo8tx+Y+dou9zQBBhnADADVg6IPzLesLpgyVy+WyuRsguBBuAMDL3l212bJ+z5DLFR5W2+ZugOBDuAEAL5u3dJ1brVXzhurR6TQHugGCD+EGALyouLjEsj5t3A02dwIEL8INAHhR6j2z3Wp339zTgU6A4EW4AQAv2fDNT5b1i89va3MnQHAj3ACAFxhj9NhL71uOcXUUYC/CDQB4wbOvf2RZX/jUMJs7AUC4AYBTZIzRf7780a1+1cVnK6x2LQc6AoIb4QYATtHN4+ZY1odef6HNnQCQCDcAcEqMMTpeUORWf2z0dfY3A0AS4QYATsk7H1nfjbhd66Y2dwLgBMINAJykY8cL9dq7n7nV33hyqAPdADiBcAMAJ8nTuTYR4Tw/CnAS4QYAToIxxrJ+2w0X2dwJgD8i3ADASeg7epZlPeXCs2zuBMAfEW4AoJreXWV9EvH8x4fY3AkAK4QbAKiGA4ePat7SdZZjdSPDbe4GgBXCDQBUw4iHXresL5k+wuZOAHhCuAGAKjp2vNCyPm5YLx6OCfgQwg0AVNH905e61dq2aqJu5yTa3gsAzxwPNzNmzFBiYqIiIiLUvXt3rV+/vsLpjxw5opEjR6pZs2YKDw/XGWecoeXLl9vULYBg9cu+I9q195BbPX3MXxzoBkBFHH1c7eLFizV27FjNnDlT3bt31/Tp05WSkqItW7aoSZMmbtMXFhbq8ssvV5MmTbRkyRLFx8frp59+UmxsrP3NAwgqdz+6yLLO4SjA9zgabqZNm6bhw4dryJBfL5+cOXOmli1bpjlz5mj8+PFu08+ZM0eHDh3S2rVrVbv2r3cATUxMrPA9CgoKVFBQUPY6JyfHex8AQFCwOhwlSY+PZa8N4IscOyxVWFioDRs2KDk5+bdmQkKUnJysdeusL7N89913lZSUpJEjRyouLk5nn322HnvsMZWUlHh8n/T0dMXExJT9JCQkeP2zAAhchUXF2rIj060e1zBabVvFOdARgMo4Fm4OHDigkpISxcWV3zjExcUpM9N9QyJJ27dv15IlS1RSUqLly5dr4sSJmjp1qh555BGP7zNhwgRlZ2eX/ezevdurnwNAYPtkw4+W9Rcm9be5EwBV5ehhqeoqLS1VkyZN9NJLLyk0NFRdunTRL7/8oilTpigtLc1ynvDwcIWHc2MtACfnhYUfu9VeeWSQA50AqCrHwk2jRo0UGhqqrKyscvWsrCw1bdrUcp5mzZqpdu3aCg0NLau1b99emZmZKiwsVFhYWI32DCC4HM0vsKzHRtWxuRMA1eHYYamwsDB16dJFGRkZZbXS0lJlZGQoKSnJcp4LLrhAW7duVWlpaVnthx9+ULNmzQg2ALzu4y9+cKslxjdyoBMA1eHofW7Gjh2r2bNna968efruu+90xx13KC8vr+zqqUGDBmnChAll099xxx06dOiQRo0apR9++EHLli3TY489ppEjRzr1EQAEsO+2u5//9+Q9f3WgEwDV4eg5N6mpqdq/f78mTZqkzMxMderUSStWrCg7yXjXrl0KCfktfyUkJOiDDz7QmDFjdO655yo+Pl6jRo3SuHHjnPoIAAKUMUbrNm1zq4eGOn7vUwCVcBljjNNN2CknJ0cxMTHKzs5WdHS00+0A8FEjJy9Q5oHy98Vq2ihaMyZylRTghOr8/ea/IADwB6WlpW7BRpIu79HBgW4AVBfhBgD+4IYxL1nW+/TsaHMnAE4G4QYAfifvmPXl3y+mDeA5UoCfINwAwO/c99Q/3WrhYbXVpEGUA90AOBmEGwD4P7l5xy3PtXnt8SEOdAPgZBFuAEC/nkR8y/1zLce4/BvwL/zGAoCk/n9/xbL+YtoAmzsBcKoINwCCXt6xAhUVl1iOca4N4H8INwCC3pIPNlrWH7z9aps7AeANjj5+AQCclpt3XO+u2uxWnzbuBrVq3tCBjgCcKvbcAAhqnk4iJtgA/otwAyBovfOR+x4bSbrjpkts7gSANxFuAASt+e+ss6wnJ7W3uRMA3kS4ARCUPtu83bL+z2dut7kTAN5GuAEQlKbM+bdb7aarujrQCQBvq1a4GTRokHJzc8teb968WUVFRV5vCgBq0hf/22lZ73XhWfY2AqBGVCvcvPHGGzp27FjZ64suuki7d+/2elMAUJMen73CrXbNJecqqm6EA90A8LZq3efGGFPhawDwZbv2HtI/XnjPcmxA7242dwOgpnATPwBBYetP+zRu2r8sx1KvPF9htdkcAoGi2r/N3377rTIzMyX9uufm+++/19GjR8tNc+6553qnOwDwgn2Hcj0GG0m6sdf5NnYDoKZVO9xcdtll5Q5HXXPNNZIkl8slY4xcLpdKSqwfQAcAdisoLNIdD7/hcXzy3X1s7AaAHaoVbnbs2FFTfQBAjej/91c8jj0+9i9q2yrOxm4A2KFa4aZVq1Y11QcAeN2h7DyPY9ysDwhcJ3UG3Y8//qh33nlHO3fulMvlUuvWrXXdddepTZs23u4PAE6KMUbDJ71mOUawAQJbtcNNenq6Jk2apNLSUjVp0kTGGO3fv1/jx4/XY489pnvvvbcm+gSAKjPGqO/oWZZjS6aPsLkbAHar1k38Vq1apQcffFAPPPCADhw4oL179yozM7Ms3IwfP17/+c9/aqpXAKiSZ1//yLIe1zBaLpfL5m4A2M1lqnEnvtTUVMXGxmrWLOv/Ed12223Kzc3VwoULvdagt+Xk5CgmJkbZ2dmKjo52uh0AXlZcXKLUe2Zbjr319G0KCeGReoA/qs7f72r9lq9fv14DBw70OD5w4EB99tln1VkkAHjV5//daVl/7fFbCTZAkKjWb3pWVpYSExM9jrdu3brsBn8A4ISXl3zqXps8SHUiwxzoBoATqhVujh8/rrAwzxuI2rVrq7Cw8JSbAoCTlXP0mFutfnQdBzoB4JRqXy318ssvq169epZjubm5p9wQAJysH3ZmOd0CAB9QrXDTsmVLzZ5tfaLe76cBALsZYzTh6bfd6k/9va8D3QBwUrXCzc6dO2uoDQA4Nfc8ucSy3rpFI5s7AeC0ap1z89FHH6lDhw7KyclxG8vOztZZZ52lTz75xGvNAUBVZHz2nX7ac9CtPvDaPznQDQCnVSvcTJ8+XcOHD7e8vjwmJkYjRozQtGnTvNYcAFTFCws/tqxfd1knexsB4BOqFW42b96sXr16eRy/4oortGHDhlNuCgBO1ZP3XO90CwAcUu373NSuXdvjeK1atbR///5TbgoAqupwTr5b7c/d2um0lo0d6AaAL6hWuImPj9f//vc/j+Nff/21mjVrdspNAUBVjXn8Tbfa0L9e4EAnAHxFtcLNVVddpYkTJ+r48eNuY8eOHVNaWpquueYarzUHABX5ac8h5ea5b4+4GzEQ3Kr14MysrCx17txZoaGhuuuuu9SuXTtJ0vfff68ZM2aopKREGzduVFxcXI01fKp4cCYQGL7dtlcTn33Hre6StOSZ2+1vCECNqs7f72rd5yYuLk5r167VHXfcoQkTJuhELnK5XEpJSdGMGTN8OtgACAwrPvlGs5dY33bilUcG29wNAF9T7ccvtGrVSsuXL9fhw4e1detWGWPUtm1b1a9fvyb6A4ByiotLPAabCzqfrpioSJs7AuBrqh1uTqhfv766du3qzV4AoFKz3rQONue1T9CYQZfZ3A0AX3TS4QYA7FZUVKKPPv/erd77z+fqlr/0cKAjAL6oWldLAYCTbrrX+sG9BBsAv0e4AeAXCouKLeuT7+5jcycAfB3hBoBfePXttZb1Dqdx41AA5RFuAPiFf6/51q22ZPoIBzoB4OsINwD8lsvlcroFAD6IcAPA51ndSL154xgHOgHgDwg3AHxaYVGx+o6e5VYfkXqxA90A8AeEGwA+a9feQ+p378uWY21bNbG5GwD+gnADwCcdO16oMY+/6XE8PKy2jd0A8CeEGwA+p7S0VDePm+Nx/I0nh9rYDQB/Q7gB4HPueXKJx7GZaQMUEc5eGwCe8WwpAD6loLBIu/Yeshz75zO329wNAH/kE3tuZsyYocTEREVERKh79+5av359leZbtGiRXC6XrrvuupptEIBtPN2JmBv2Aagqx8PN4sWLNXbsWKWlpWnjxo3q2LGjUlJStG/fvgrn27lzp+69915ddNFFNnUKwA4frv3OrTbnkcHcsA9AlTkebqZNm6bhw4dryJAh6tChg2bOnKk6depozhzPJxOWlJRowIABevjhh9WmTRsbuwVQk4qLSyzrMVGRNncCwJ85Gm4KCwu1YcMGJScnl9VCQkKUnJysdevWeZzvH//4h5o0aaKhQyu/YqKgoEA5OTnlfgD4ptR7ZrvV/nJZJ/sbAeDXHA03Bw4cUElJieLi4srV4+LilJmZaTnPp59+qldeeUWzZ7tvBK2kp6crJiam7CchIeGU+wbgfXv2HbGsX39FZ3sbAeD3HD8sVR25ubkaOHCgZs+erUaNGlVpngkTJig7O7vsZ/fu3TXcJYCT8bdHF1nWIyPCbO4EgL9z9FLwRo0aKTQ0VFlZWeXqWVlZatq0qdv027Zt086dO9W7d++yWmlpqSSpVq1a2rJli0477bRy84SHhys8PLwGugfgLbszD1vWF08dbnMnAAKBo3tuwsLC1KVLF2VkZJTVSktLlZGRoaSkJLfpzzzzTP33v//Vpk2byn6uvfZaXXrppdq0aROHnAA/NTp9sVstOam9atUKdaAbAP7O8Zv4jR07VoMHD9b555+vbt26afr06crLy9OQIUMkSYMGDVJ8fLzS09MVERGhs88+u9z8sbGxkuRWB+Afvt2217J+x02X2NwJgEDheLhJTU3V/v37NWnSJGVmZqpTp05asWJF2UnGu3btUkiIX50aBKAaJj77jlvt3DNaONAJgEDhMsYYp5uwU05OjmJiYpSdna3o6Gin2wGC2v9b9bXmLnW/I/GS6SO4aR+Acqrz95tdIgAcsezj/1oGm84dWhJsAJwSwg0A2xUWFWvOv9ZYjo0f1svmbgAEGsINANuNnLzQsn795Z0VGspmCcCpcfyEYgDB51B2nlvtrv6X6tLu7RzoBkCg4b9IAGxlFWwkEWwAeA3hBoCtnp630q1239AUBzoBEKgINwBsZXXTvu7ntnagEwCBinADAAACCuEGgG2uHzXTrXZ5j/YOdAIgkBFuANgi71iBZf2mq7ra3AmAQEe4AVDjiotLNGj8q5ZjsVF1bO4GQKAj3ACoUcYYpd4z23LsjSeH2twNgGBAuAFQo/qOnmVZ735ua0WE17a5GwDBgHADoMZs+OYnj2Pc2wZATSHcAKgRxhg99tL7lmNzH7vF3mYABBXCDYAa8cLCjy3rzz1wk6LqRtjcDYBgQrgBUCM++vx7t1rP7meqeZNY+5sBEFQINwC8LufoMcv6nf0usbkTAMGIcAPA625Le92tNuLGi+VyuRzoBkCwIdwA8KojufkqKi5xq/fs3s6BbgAEI8INAK8xxmjog/Pd6g1j66pWrVAHOgIQjAg3ALxmyb83WtafHn+jzZ0ACGaEGwBes2j5F261Lh1aqW5kuAPdAAhWhBsAXnEkN9+yfv+IK23uBECwI9wA8Aqrc22G973IgU4ABDvCDYBTNvrxNy3rvS46y+ZOAIBwA+AUHc0v0O69h9zqDWPrOtANABBuAJyiwRNetazPeuhmmzsBgF8RbgCctAOHj1rWn3vgJu5GDMAxhBsAJ23EQ+6PWZDEwzEBOIpwA8CrFk8d7nQLAIIc4QbASdl/KNet1q51Ux6zAMBxhBsAJ+V/P+5xqz3ADfsA+ADCDYCT8vyCVW41HrMAwBcQbgBUm9VVUo3rRznQCQC4I9wAqDarq6Rat2joQCcA4I5wA6BaiotLLOv3DU2xuRMAsEa4AVAtk57/f261i89vy037APgMwg2AKlvz1TZt2ZHpVr/75p4OdAMA1gg3AKok80COps390HKMvTYAfAnhBkCl8o4VaOTkBZZjk+/uY3M3AFAxwg2ASt364HzLes/uZ6rDac1s7gYAKka4AVChzAM5lldIxUbV0Z39LnGgIwCoWC2nGwDgm4wxmvziMm3e8rPl+CuPDLK5IwCoGvbcALA0/53PPAabVx8dbHM3AFB1hBsAlt5dtdmyfu4ZLRRdL9LmbgCg6jgsBcBNbt5xy/oVF3TQiBsvtrkbAKgewg0AN7fcP9et1u/qbup7RWf7mwGAauKwFIByhk96zbJ+Xc+ONncCACeHcAOgzD8/3KhD2XmWY7VqhdrcDQCcHMINAEnSoew8LXhvveXYnEe4OgqA/yDcANDR/AKPh6Nm/2OgYqK4OgqA/yDcANDgCa9a1ifdeY0axNS1uRsAODVcLQUEMWOMJj33ruVY6xaN1LFdC5s7AoBTR7gBgtjfHl2kvfuzLcem3Hu9zd0AgHdwWAoIYp6CzcKnhsnlctncDQB4h0+EmxkzZigxMVERERHq3r271q+3vmJDkmbPnq2LLrpI9evXV/369ZWcnFzh9ACszX9nnWW9b0oXhdVmpy4A/+V4uFm8eLHGjh2rtLQ0bdy4UR07dlRKSor27dtnOf3q1avVr18/rVq1SuvWrVNCQoKuuOIK/fLLLzZ3Dvi3dz5yf3bUAyOuUr+rujrQDQB4j8sYY5xsoHv37uratauef/55SVJpaakSEhL0t7/9TePHj690/pKSEtWvX1/PP/+8Bg0aVOn0OTk5iomJUXZ2tqKjo0+5f8Af3TB6lkotfvX/+cztDnQDAJWrzt9vR/fcFBYWasOGDUpOTi6rhYSEKDk5WevWWe8y/6P8/HwVFRWpQYMGluMFBQXKyckp9wMEs+LiEstgM/sfAx3oBgC8z9Fwc+DAAZWUlCguLq5cPS4uTpmZmVVaxrhx49S8efNyAen30tPTFRMTU/aTkJBwyn0D/irvWIFS75ltOcb9bAAECsfPuTkVjz/+uBYtWqS3335bERERltNMmDBB2dnZZT+7d++2uUvANxQWFWvQeOub9b0wqb/N3QBAzXH0kohGjRopNDRUWVlZ5epZWVlq2rRphfM+9dRTevzxx7Vy5Uqde+65HqcLDw9XeHi4V/oF/Fm/e1/2OBbXkPPPAAQOR/fchIWFqUuXLsrIyCirlZaWKiMjQ0lJSR7ne/LJJzV58mStWLFC559/vh2tAn7t/ulLPY4temq4fY0AgA0cv5nF2LFjNXjwYJ1//vnq1q2bpk+frry8PA0ZMkSSNGjQIMXHxys9PV2S9MQTT2jSpElasGCBEhMTy87NqVevnurVq+fY5wB81ZqvtmnLDutz2JZMH8HN+gAEHMfDTWpqqvbv369JkyYpMzNTnTp10ooVK8pOMt61a5dCQn7bwfTiiy+qsLBQffv2LbectLQ0PfTQQ3a2Dvi891Z/rVffXms5RrABEKgcv8+N3bjPDYLF6+9+prczNlmOzf7HQK6OAuBX/OY+NwBqxpqvtnkMNnff3JNgAyCgEW6AALPj5wOaNvdDy7GLurTVJV3PsLkjALCX4+fcAPCO0tJS3TDmJY/jd9/ck2ADICiw5wYIAAcOH60w2NxyXQ+CDYCgQbgB/NyBw0c14qHXPY4nNGug3pd6vtElAAQaDksBfqy4uKTCYHPTVV11Q0oXGzsCAOcRbgA/tnbTNo9j3McGQLAi3AB+7JnXPrKs//OZ223uBAB8B+fcAH7oeEGRrh8103Js8VSeFQUguBFuAD9z7HihBtz3iuVYygVnqVatUJs7AgDfQrgB/IgxRjePm+NxfPgNF9rYDQD4JsIN4Ef6jp7lcWzhU8M4gRgAxAnFgN8YP+1fHse4MgoAfkO4AXzcrr2HNG3uh9qdedhynCujAKA8wg3gw4wxGvP4mx7Hl0wfYWM3AOAfOOcG8FHGmArPsXkxbQCHogDAAuEG8EHfbttbYbAZ1vdCNWkQZWNHAOA/OCwF+BhjjCY++47H8TmPDFZMVKSNHQGAf2HPDeBDKjsU9cojgwg2AFAJwg3gQ4Y8MM/j2NT7+io2qo6N3QCAf+KwFOADjDH626OLlJt33HKcy70BoOoIN4DDKjsU9foTt9rYDQD4Pw5LAQ4qLi6pMNi89PDNiowIs7EjAPB/hBvAIUfzC5R6z2yP4w/f1VsNY+vZ2BEABAYOSwEOyDyQo5GTF3gcn/XQzWpUn2ADACeDcAPY7MefsjR+2tsex2f/Y6AaxNS1sSMACCyEG8BGlQWbt56+TSEhHC0GgFNBuAFsYIzREy9/oC/+t9PjNNPG3UCwAQAvINwANayyS72bNY7R8w/2s7EjAAhs/DcRqGEVBZv4JrEEGwDwMvbcADXouTdWeRy76uKzNfT6C23sBgCCA+EGqAHfb8/UA88s9Tg+bdwNatW8oX0NAUAQIdwAXpZ/rLDCYPPmtNsUGsoRYQCoKYQbwAveXPGl3v/kG+UcPVbhdOOG9SLYAEANI9wAp+jGsS+ppKS00ummT0hVQtP6NnQEAMGNcAOcgoXL1lcp2Pzzmdtt6AYAIBFugJN2/aiZlU4T1zBaMyZyqTcA2IlwA1RTbt5x3XL/XI/j11/eWZ07tFSzxjGKiYq0rzEAgCTCDVAte/dn665HFnocv/WvF+jqS86xsSMAwB8RboBqqCjYPDH2rzq9VRMbuwEAWCHcAJUwxujtlZv0xnufe5zm2QduUnyTWPuaAgB4RLgBKnAoO0/DJ71W4TSzHrpZjerXs6kjAEBlCDeAB8++/pE+/uKHCqdZ9NRw1a4dalNHAICqINwAf/D/Vn2tuUvXVjrd1Pv6EmwAwAcRboDf2fDNT5UGm6svOUd9enZUw1gORQGALyLcAPr1pOHJLy7T5i0/VzjdC5P6K65htE1dAQBOBuEGQc8Yo76jZ1U4zYgbL9YVF3SwqSMAwKkg3CAoHC8o0s5fDqqktPxzoDZ+u0tLMzZ5nK/Heafp9tSLVTcyvIY7BAB4C+EGAWvfoVwdOHxUE59956Tmv+riszX0+gu93BUAoKYRbhCQ+o6aKXMK889MG6DGDaK81g8AwD6EG/itwqJi7dpzSD/tPahvt2Vq9fotXlnuG08OVUR4ba8sCwBgP8IN/E5xcYkGjn9VhUXFJzV/aGiIIn8XXo7mF0j69aGXV118tlwul1f6BAA4g3ADn2CM0S/7jqigoHxg+eGnLG3bvV+rPt+iOhFhyj9eeErvs3jqcNWqxY33ACCQEW5Qo4wxWv6f/ynrYI7C/i9UFBaXaNnH/1V0vUhFhtfW0fwC5R0rqHRZJxtsXkwboCacPwMAQYNwg2rbsiNT6zZtV1htz1+f0tJSvV3BJdaSlHP0mHKOHvNqb2e3ba7mTWLV7ZzWio2KVGJ8Qw4zAUCQIdygzP5DuVr9xQ8qKCjSus3blXkgR80bx5SbZs/+bIe682xe+hDVq8N9aAAAv/KJcDNjxgxNmTJFmZmZ6tixo5577jl169bN4/RvvfWWJk6cqJ07d6pt27Z64okndNVVV9nYsX/4dMNW7fzlgI4VFGnFp9+oaaNohYaEuE1XVFyifYdyLZfhVJip/bvzYoqKSyT9ulembcsmOvO0ZooMr622rZpUuPcIABCcHP/LsHjxYo0dO1YzZ85U9+7dNX36dKWkpGjLli1q0qSJ2/Rr165Vv379lJ6ermuuuUYLFizQddddp40bN+rss8924BP8Kjv3mD7bvF37D+Xq7YxNalS/niLCnLuc+Oesw261zAM5DnRS3oVdTlfI/x0mOpJzTK1bNFT705pJkmqFhqptqybshQEAnBKXMeZU7nV2yrp3766uXbvq+eefl/TruRoJCQn629/+pvHjx7tNn5qaqry8PL333ntltT/96U/q1KmTZs6cWen75eTkKCYmRtnZ2YqO9s4DEIuLSzTm8Td98pBNTevYroWi6kV4HM89elzxcbG66aquPMIAAHDSqvP329E9N4WFhdqwYYMmTJhQVgsJCVFycrLWrVtnOc+6des0duzYcrWUlBQtXbrUcvqCggIVFPx2JU5Ojvf3Xmz/+UBABps/ndtaOXnH1euisxUaUv6k3LqR4WrXOo7DQgAAn+PoX6YDBw6opKREcXFx5epxcXH6/vvvLefJzMy0nD4zM9Ny+vT0dD388MPeadiD4wVFNbp8b0jqdJoOHM5VxzMTlNi8oeU0ISEund6ysRrG1rO5OwAAvCfg/9s9YcKEcnt6cnJylJCQUOPve0ZinFrHN9JZbZvX+HtVpEVcrFo2a8Dl0ACAoOFouGnUqJFCQ0OVlZVVrp6VlaWmTZtaztO0adNqTR8eHq7w8Jo91+OcM+K1ZPqIcjXCBAAAznC/LthGYWFh6tKlizIyMspqpaWlysjIUFJSkuU8SUlJ5aaXpA8//NDj9HZwuVxuPwAAwBmOH5YaO3asBg8erPPPP1/dunXT9OnTlZeXpyFDhkiSBg0apPj4eKWnp0uSRo0apUsuuURTp07V1VdfrUWLFunLL7/USy+95OTHAAAAPsLxcJOamqr9+/dr0qRJyszMVKdOnbRixYqyk4Z37dqlkN/deK5Hjx5asGCBHnzwQd1///1q27atli5d6ug9bgAAgO9w/D43dquJ+9wAAICaVZ2/346ecwMAAOBthBsAABBQCDcAACCgEG4AAEBAIdwAAICAQrgBAAABhXADAAACCuEGAAAEFMINAAAIKI4/fsFuJ27InJOT43AnAACgqk783a7KgxWCLtzk5uZKkhISEhzuBAAAVFdubq5iYmIqnCboni1VWlqqPXv2KCoqSi6Xy+l2qi0nJ0cJCQnavXs3z8YS6+P3WBflsT7KY338hnVRnr+sD2OMcnNz1bx583IP1LYSdHtuQkJC1KJFC6fbOGXR0dE+/SW0G+vjN6yL8lgf5bE+fsO6KM8f1kdle2xO4IRiAAAQUAg3AAAgoBBu/Ex4eLjS0tIUHh7udCs+gfXxG9ZFeayP8lgfv2FdlBeI6yPoTigGAACBjT03AAAgoBBuAABAQCHcAACAgEK4AQAAAYVw46CSkhJNnDhRrVu3VmRkpE477TRNnjy57LkZRUVFGjdunM455xzVrVtXzZs316BBg7Rnz54Kl/vQQw/J5XKV+znzzDPt+EgnrbJ1IUm33HKL2+fq1atXpcueMWOGEhMTFRERoe7du2v9+vU1+VG8oirr44/r4sTPlClTPC7XH78bJ+Tm5mr06NFq1aqVIiMj1aNHD33xxRdl48YYTZo0Sc2aNVNkZKSSk5P1448/Vrpcf/x+VLQugmm7cUJl341g2nZUti6CZrth4JhHH33UNGzY0Lz33ntmx44d5q233jL16tUzzzzzjDHGmCNHjpjk5GSzePFi8/3335t169aZbt26mS5dulS43LS0NHPWWWeZvXv3lv3s37/fjo900ipbF8YYM3jwYNOrV69yn+vQoUMVLnfRokUmLCzMzJkzx3zzzTdm+PDhJjY21mRlZdX0RzolVVkfv18Pe/fuNXPmzDEul8ts27bN43L98btxwo033mg6dOhgPv74Y/Pjjz+atLQ0Ex0dbX7++WdjjDGPP/64iYmJMUuXLjWbN2821157rWndurU5duyYx2X66/ejonURTNuNEyr7bgTTtqOydREs2w3CjYOuvvpqc+utt5ar/fWvfzUDBgzwOM/69euNJPPTTz95nCYtLc107NjRW23aoirrYvDgwaZPnz7VWm63bt3MyJEjy16XlJSY5s2bm/T09FPqt6adzHejT58+pmfPnhUu1x+/G8YYk5+fb0JDQ817771Xrt65c2fzwAMPmNLSUtO0aVMzZcqUsrEjR46Y8PBws3DhQo/L9cfvR2XrwkqgbjeMqdr6CJZtx8l8NwJ1u8FhKQf16NFDGRkZ+uGHHyRJmzdv1qeffqorr7zS4zzZ2dlyuVyKjY2tcNk//vijmjdvrjZt2mjAgAHatWuXN1v3uqqui9WrV6tJkyZq166d7rjjDh08eNDjMgsLC7VhwwYlJyeX1UJCQpScnKx169bVzAfxkup+N7KysrRs2TINHTq00mX723dDkoqLi1VSUqKIiIhy9cjISH366afasWOHMjMzy/1bx8TEqHv37h7/rf31+1HZurASqNsNqerrIxi2HdX9bgT0dsPpdBXMSkpKzLhx44zL5TK1atUyLpfLPPbYYx6nP3bsmOncubPp379/hctdvny5efPNN83mzZvNihUrTFJSkmnZsqXJycnx9kfwmqqsi4ULF5p33nnHfP311+btt9827du3N127djXFxcWWy/zll1+MJLN27dpy9b///e+mW7duNfZZvKG6340nnnjC1K9fv8JDMMb453fjhKSkJHPJJZeYX375xRQXF5vXXnvNhISEmDPOOMOsWbPGSDJ79uwpN88NN9xgbrzxRsvl+fP3o6J18UeBvN04obL1EUzbjup8NwJ5u0G4cdDChQtNixYtzMKFC83XX39t5s+fbxo0aGDmzp3rNm1hYaHp3bu3Oe+880x2dna13ufw4cMmOjravPzyy95q3euqsy5O2LZtm5FkVq5caTnuzxuo6q6Pdu3ambvuuqva7+MP340Ttm7dai6++GIjyYSGhpquXbuaAQMGmDPPPDPowk1F6+L3An27cUJV18cJgbztqM66COTtBoelHPT3v/9d48eP10033aRzzjlHAwcO1JgxY5Senl5uuqKiIt1444366aef9OGHH1b7kfSxsbE644wztHXrVm+271VVXRe/16ZNGzVq1Mjj52rUqJFCQ0OVlZVVrp6VlaWmTZt6tX9vq876+OSTT7RlyxYNGzas2u/jD9+NE0477TR9/PHHOnr0qHbv3q3169erqKhIbdq0Kfv3rM6/tT9/PypaFycEw3bjhKqsj98L5G1HVddFoG83CDcOys/PV0hI+X+C0NBQlZaWlr0+sYH68ccftXLlSjVs2LDa73P06FFt27ZNzZo1O+Wea0pV1sUf/fzzzzp48KDHzxUWFqYuXbooIyOjrFZaWqqMjAwlJSV5p/EaUp318corr6hLly7q2LFjtd/HH74bf1S3bl01a9ZMhw8f1gcffKA+ffqodevWatq0abl/65ycHH3++ece/639+ftxgtW6kIJnu/FHntbHHwXytuOEytZFwG83nN51FMwGDx5s4uPjyy73/de//mUaNWpk7rvvPmPMr7uUr732WtOiRQuzadOmcpfhFRQUlC2nZ8+e5rnnnit7fc8995jVq1ebHTt2mDVr1pjk5GTTqFEjs2/fPts/Y1VVti5yc3PNvffea9atW2d27NhhVq5caTp37mzatm1rjh8/XracP66LRYsWmfDwcDN37lzz7bffmttuu83ExsaazMxM2z9jdVS2Pk7Izs42derUMS+++KLlcgLhu3HCihUrzPvvv2+2b99u/v3vf5uOHTua7t27m8LCQmPMr5eCx8bGlp1b0adPH7dLwQPl+1HRugim7cYJFa2PYNt2VPZ7YkxwbDcINw7Kyckxo0aNMi1btjQRERGmTZs25oEHHijbAO3YscNIsvxZtWpV2XJatWpl0tLSyl6npqaaZs2ambCwMBMfH29SU1PN1q1bbf501VPZusjPzzdXXHGFady4saldu7Zp1aqVGT58uNuG5o/rwhhjnnvuOdOyZUsTFhZmunXrZj777DO7PtZJq2x9nDBr1iwTGRlpjhw5YrmcQPhunLB48WLTpk0bExYWZpo2bWpGjhxZ7nOXlpaaiRMnmri4OBMeHm4uu+wys2XLlnLLCJTvR0XrIpi2GydUtD6CbdtR2e+JMcGx3XAZ87tbngIAAPg5zrkBAAABhXADAAACCuEGAAAEFMINAAAIKIQbAAAQUAg3AAAgoBBuAABAQCHcAACAgEK4AeB3XC6Xli5d6nQbAHwU4QaAm/379+uOO+5Qy5YtFR4erqZNmyolJUVr1qxxujUAqFQtpxsA4Huuv/56FRYWat68eWrTpo2ysrKUkZGhgwcPOt0aAFSKPTcAyjly5Ig++eQTPfHEE7r00kvVqlUrdevWTRMmTNC1115bNt20adN0zjnnqG7dukpISNCdd96po0ePlo3PnTtXsbGxeu+999SuXTvVqVNHffv2VX5+vubNm6fExETVr19fd999t0pKSsrmS0xM1OTJk9WvXz/VrVtX8fHxmjFjRoU97969WzfeeKNiY2PVoEED9enTRzt37vQ4/erVq+VyufTBBx/ovPPOU2RkpHr27Kl9+/bp/fffV/v27RUdHa3+/fsrPz+/bL7S0lKlp6erdevWioyMVMeOHbVkyZKy8ZKSEg0dOrRsvF27dnrmmWfKvfctt9yi6667Tk899ZSaNWumhg0bauTIkSoqKqr03wZA1RBuAJRTr1491atXT0uXLlVBQYHH6UJCQvTss8/qm2++0bx58/TRRx/pvvvuKzdNfn6+nn32WS1atEgrVqzQ6tWr9Ze//EXLly/X8uXL9dprr2nWrFnlAoIkTZkyRR07dtRXX32l8ePHa9SoUfrwww8t+ygqKlJKSoqioqL0ySefaM2aNapXr5569eqlwsLCCj/rQw89pOeff15r164tC0jTp0/XggULtGzZMv373//Wc889VzZ9enq65s+fr5kzZ+qbb77RmDFjdPPNN+vjjz+W9Gv4adGihd566y19++23mjRpku6//369+eab5d531apV2rZtm1atWqV58+Zp7ty5mjt3boW9AqgGpx9LDsD3LFmyxNSvX99ERESYHj16mAkTJpjNmzdXOM9bb71lGjZsWPb61VdfNZLM1q1by2ojRowwderUMbm5uWW1lJQUM2LEiLLXrVq1Mr169Sq37NTUVHPllVeWvZZk3n77bWOMMa+99ppp166dKS0tLRsvKCgwkZGR5oMPPrDsddWqVUaSWblyZVktPT3dSDLbtm0r129KSooxxpjjx4+bOnXqmLVr15Zb1tChQ02/fv08rpeRI0ea66+/vuz14MGDTatWrUxxcXFZ7YYbbjCpqakelwGgethzA8DN9ddfrz179ujdd99Vr169tHr1anXu3Lnc3oWVK1fqsssuU3x8vKKiojRw4EAdPHiw3GGcOnXq6LTTTit7HRcXp8TERNWrV69cbd++feXePykpye31d999Z9nr5s2btXXrVkVFRZXtdWrQoIGOHz+ubdu2Vfg5zz333HJ91KlTR23atLHsbevWrcrPz9fll19e9j716tXT/Pnzy73PjBkz1KVLFzVu3Fj16tXTSy+9pF27dpV737POOkuhoaFlr5s1a+a2DgCcPE4oBmApIiJCl19+uS6//HJNnDhRw4YNU1pamm655Rbt3LlT11xzje644w49+uijatCggT799FMNHTpUhYWFqlOnjiSpdu3a5Zbpcrksa6WlpSfd59GjR9WlSxe98cYbbmONGzeucN7f91JZbyfOJ1q2bJni4+PLTRceHi5JWrRoke69915NnTpVSUlJioqK0pQpU/T55597fN8/vg+AU0e4AVAlHTp0KLu3zIYNG1RaWqqpU6cqJOTXHcB/PK/kVHz22Wdur9u3b285befOnbV48WI1adJE0dHRXuvhjzp06KDw8HDt2rVLl1xyieU0a9asUY8ePXTnnXeW1SrbewTA+zgsBaCcgwcPqmfPnnr99df19ddfa8eOHXrrrbf05JNPqk+fPpKk008/XUVFRXruuee0fft2vfbaa5o5c6bXelizZo2efPJJ/fDDD5oxY4beeustjRo1ynLaAQMGqFGjRurTp48++eQT7dixQ6tXr9bdd9+tn3/+2Ws9RUVF6d5779WYMWM0b948bdu2TRs3btRzzz2nefPmSZLatm2rL7/8Uh988IF++OEHTZw4UV988YXXegBQNey5AVBOvXr11L17dz399NPatm2bioqKlJCQoOHDh+v++++XJHXs2FHTpk3TE088oQkTJujiiy9Wenq6Bg0a5JUe7rnnHn355Zd6+OGHFR0drWnTpiklJcVy2jp16ug///mPxo0bp7/+9a/Kzc1VfHy8LrvsMq/vyZk8ebIaN26s9PR0bd++XbGxsercuXPZehkxYoS++uorpaamyuVyqV+/frrzzjv1/vvve7UPABVzGWOM000AwAmJiYkaPXq0Ro8e7XQrAPwUh6UAAEBAIdwAAICAwmEpAAAQUNhzAwAAAgrhBgAABBTCDQAACCiEGwAAEFAINwAAIKAQbgAAQEAh3AAAgIBCuAEAAAHl/wN4YJiRmXzNlQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xbars = SimulateSample()\n",
    "cdf = thinkstats2.Cdf(xbars)\n",
    "thinkplot.Cdf(cdf)\n",
    "thinkplot.Config(xlabel='Sample mean',\n",
    "                 ylabel='CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d0a0dc3",
   "metadata": {},
   "source": [
    "In each iteration, we choose `n` values from a normal distribution with\n",
    "the given parameters, and compute the sample mean, `xbar`. We run 1000\n",
    "simulations and then compute the distribution, `cdf`, of the estimates.\n",
    "The result is shown in\n",
    "Figure [\\[estimation1\\]](#estimation1){reference-type=\"ref\"\n",
    "reference=\"estimation1\"}. This distribution is called the **sampling\n",
    "distribution** of the estimator. It shows how much the estimates would\n",
    "vary if we ran the experiment over and over.\n",
    "\n",
    "The mean of the sampling distribution is pretty close to the\n",
    "hypothetical value of $\\mu$, which means that the experiment yields the\n",
    "right answer, on average. After 1000 tries, the lowest result is 82 kg,\n",
    "and the highest is 98 kg. This range suggests that the estimate might be\n",
    "off by as much as 8 kg."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "e4c223b2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "90.06484084475981"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(xbars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "931496f2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(85.96712753834224, 94.13764998760222)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ci = cdf.Percentile(5), cdf.Percentile(95)\n",
    "ci"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "4299ce1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5378143765427614"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stderr = RMSE(xbars, 90)\n",
    "stderr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c103ec8f",
   "metadata": {},
   "source": [
    "There are two common ways to summarize the sampling distribution:\n",
    "\n",
    "-   **Standard error** (SE) is a measure of how far we expect the\n",
    "    estimate to be off, on average. For each simulated experiment, we\n",
    "    compute the error, $\\bar{x} - \\mu$, and then compute the root mean\n",
    "    squared error (RMSE). In this example, it is roughly 2.5 kg.\n",
    "\n",
    "-   A **confidence interval** (CI) is a range that includes a given\n",
    "    fraction of the sampling distribution. For example, the 90%\n",
    "    confidence interval is the range from the 5th to the 95th\n",
    "    percentile. In this example, the 90% CI is $(86, 94)$ kg."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7deaf510",
   "metadata": {},
   "source": [
    "Standard errors and confidence intervals are the source of much\n",
    "confusion:\n",
    "\n",
    "-   People often confuse standard error and standard deviation. Remember\n",
    "    that standard deviation describes variability in a measured\n",
    "    quantity; in this example, the standard deviation of gorilla weight\n",
    "    is 7.5 kg. Standard error describes variability in an estimate. In\n",
    "    this example, the standard error of the mean, based on a sample of 9\n",
    "    measurements, is 2.5 kg."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64cf014b",
   "metadata": {},
   "source": [
    "One way to remember the difference is that, as sample size\n",
    "increases, standard error gets smaller; standard deviation does not.\n",
    "\n",
    "-   People often think that there is a 90% probability that the actual\n",
    "    parameter, $\\mu$, falls in the 90% confidence interval. Sadly, that\n",
    "    is not true. If you want to make a claim like that, you have to use\n",
    "    Bayesian methods (see my book, *Think Bayes*)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9812c4f0",
   "metadata": {},
   "source": [
    "The sampling distribution answers a different question: it gives you\n",
    "a sense of how reliable an estimate is by telling you how much it\n",
    "would vary if you ran the experiment again."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7b6560b",
   "metadata": {},
   "source": [
    "It is important to remember that confidence intervals and standard\n",
    "errors only quantify sampling error; that is, error due to measuring\n",
    "only part of the population. The sampling distribution does not account\n",
    "for other sources of error, notably sampling bias and measurement error,\n",
    "which are the topics of the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9185ef3",
   "metadata": {},
   "source": [
    "## Sampling bias\n",
    "\n",
    "Suppose that instead of the weight of gorillas in a nature preserve, you\n",
    "want to know the average weight of women in the city where you live. It\n",
    "is unlikely that you would be allowed to choose a representative sample\n",
    "of women and weigh them.\n",
    "\n",
    "A simple alternative would be \"telephone sampling;\" that is, you could\n",
    "choose random numbers from the phone book, call and ask to speak to an\n",
    "adult woman, and ask how much she weighs.\n",
    "\n",
    "Telephone sampling has obvious limitations. For example, the sample is\n",
    "limited to people whose telephone numbers are listed, so it eliminates\n",
    "people without phones (who might be poorer than average) and people with\n",
    "unlisted numbers (who might be richer). Also, if you call home\n",
    "telephones during the day, you are less likely to sample people with\n",
    "jobs. And if you only sample the person who answers the phone, you are\n",
    "less likely to sample people who share a phone line."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29713752",
   "metadata": {},
   "source": [
    "If factors like income, employment, and household size are related to\n",
    "weight---and it is plausible that they are---the results of your survey\n",
    "would be affected one way or another. This problem is called **sampling\n",
    "bias** because it is a property of the sampling process.\n",
    "\n",
    "This sampling process is also vulnerable to self-selection, which is a\n",
    "kind of sampling bias. Some people will refuse to answer the question,\n",
    "and if the tendency to refuse is related to weight, that would affect\n",
    "the results.\n",
    "\n",
    "Finally, if you ask people how much they weigh, rather than weighing\n",
    "them, the results might not be accurate. Even helpful respondents might\n",
    "round up or down if they are uncomfortable with their actual weight. And\n",
    "not all respondents are helpful. These inaccuracies are examples of\n",
    "**measurement error**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a97787d",
   "metadata": {},
   "source": [
    "When you report an estimated quantity, it is useful to report standard\n",
    "error, or a confidence interval, or both, in order to quantify sampling\n",
    "error. But it is also important to remember that sampling error is only\n",
    "one source of error, and often it is not the biggest."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8fee736",
   "metadata": {},
   "source": [
    "## Exponential distributions\n",
    "\n",
    "Let's play one more round of the estimation game. *I'm thinking of a\n",
    "distribution.* It's an exponential distribution, and here's a sample:\n",
    "\n",
    "`[5.384, 4.493, 19.198, 2.790, 6.122, 12.844]`\n",
    "\n",
    "What do you think is the parameter, $\\lambda$, of this distribution?\n",
    "\n",
    "In general, the mean of an exponential distribution is $1/\\lambda$, so\n",
    "working backwards, we might choose $$L= 1 / \\bar{x}$$ $L$ is an estimator\n",
    "of $\\lambda$. And not just any estimator; it is also the maximum\n",
    "likelihood estimator (see\n",
    "<http://wikipedia.org/wiki/Exponential_distribution#Maximum_likelihood>).\n",
    "So if you want to maximize your chance of guessing $\\lambda$ exactly,\n",
    "$L$ is the way to go.\n",
    "\n",
    "But we know that $\\bar{x}$ is not robust in the presence of outliers, so\n",
    "we expect $L$ to have the same problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92c023db",
   "metadata": {},
   "source": [
    "We can choose an alternative based on the sample median. The median of\n",
    "an exponential distribution is $\\ln(2) / \\lambda$, so working backwards\n",
    "again, we can define an estimator $$L_m= \\ln(2) / m$$ where $m$ is the\n",
    "sample median.\n",
    "\n",
    "To test the performance of these estimators, we can simulate the\n",
    "sampling process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9d6dbd28",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Estimate3(n=7, m=1000):\n",
    "    lam = 2\n",
    "\n",
    "    means = []\n",
    "    medians = []\n",
    "    for _ in range(m):\n",
    "        xs = np.random.exponential(1.0/lam, n)\n",
    "        L = 1 / np.mean(xs)\n",
    "        Lm = np.log(2) / thinkstats2.Median(xs)\n",
    "        means.append(L)\n",
    "        medians.append(Lm)\n",
    "\n",
    "    print('rmse L', RMSE(means, lam))\n",
    "    print('rmse Lm', RMSE(medians, lam))\n",
    "    print('mean error L', MeanError(means, lam))\n",
    "    print('mean error Lm', MeanError(medians, lam))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a550684",
   "metadata": {},
   "source": [
    "When I run this experiment with $\\lambda=2$, the RMSE of $L$ is 1.1. For\n",
    "the median-based estimator $L_m$, RMSE is 1.8. We can't tell from this\n",
    "experiment whether $L$ minimizes MSE, but at least it seems better than\n",
    "$L_m$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "604bb600",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse L 1.1170197988041468\n",
      "rmse Lm 1.7766765260102633\n",
      "mean error L 0.36424219588529966\n",
      "mean error Lm 0.5135627068422906\n"
     ]
    }
   ],
   "source": [
    "Estimate3()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a7e05a2",
   "metadata": {},
   "source": [
    "Sadly, it seems that both estimators are biased. For $L$ the mean error\n",
    "is 0.33; for $L_m$ it is 0.45. And neither converges to 0 as `m`\n",
    "increases.\n",
    "\n",
    "It turns out that $\\bar{x}$ is an unbiased estimator of the mean of the\n",
    "distribution, $1 / \\lambda$, but $L$ is not an unbiased estimator of\n",
    "$\\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a72892f",
   "metadata": {},
   "source": [
    "## Glossary\n",
    "\n",
    "-   **estimation**: The process of inferring the parameters of a\n",
    "    distribution from a sample.\n",
    "\n",
    "-   **estimator**: A statistic used to estimate a parameter.\n",
    "\n",
    "-   **mean squared error (MSE)**: A measure of estimation error.\n",
    "\n",
    "-   **root mean squared error (RMSE)**: The square root of MSE, a more\n",
    "    meaningful representation of typical error magnitude.\n",
    "\n",
    "-   **maximum likelihood estimator (MLE)**: An estimator that computes\n",
    "    the point estimate most likely to be correct.\n",
    "\n",
    "-   **bias (of an estimator)**: The tendency of an estimator to be above\n",
    "    or below the actual value of the parameter, when averaged over\n",
    "    repeated experiments.\n",
    "\n",
    "-   **sampling error**: Error in an estimate due to the limited size of\n",
    "    the sample and variation due to chance.\n",
    "\n",
    "-   **sampling bias**: Error in an estimate due to a sampling process\n",
    "    that is not representative of the population.\n",
    "\n",
    "-   **measurement error**: Error in an estimate due to inaccuracy\n",
    "    collecting or recording data.\n",
    "\n",
    "-   **sampling distribution**: The distribution of a statistic if an\n",
    "    experiment is repeated many times.\n",
    "\n",
    "-   **standard error**: The RMSE of an estimate, which quantifies\n",
    "    variability due to sampling error (but not other sources of error).\n",
    "\n",
    "-   **confidence interval**: An interval that represents the expected\n",
    "    range of an estimator if an experiment is repeated many times."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "697621aa",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "700f4df0",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose you draw a sample with size n=10 from an exponential distribution with λ=2. Simulate this experiment 1000 times and plot the sampling distribution of the estimate L. Compute the standard error of the estimate and the 90% confidence interval.\n",
    "\n",
    "Repeat the experiment with a few different values of `n` and make a plot of standard error versus `n`.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "71beaf3d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "standard error 0.8123273990929318\n",
      "confidence interval (1.2708669507738428, 3.5959240022479264)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.8123273990929318"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHHCAYAAAC2rPKaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABGNUlEQVR4nO3deVhUhf7H8c+woyJugBsK7oJbrqGZdUNNzZuVV68t4pJZqVejbrlrm2Sl2WJupdYtr5aZLeZ+NSvtai6VN7PcTQPFBRAVZDi/P/w5hWdQUJgzzLxfz8PzNN85Z+bDSPLxrDbDMAwBAAB4IB+rAwAAABQXig4AAPBYFB0AAOCxKDoAAMBjUXQAAIDHougAAACPRdEBAAAei6IDAAA8FkUHAAB4LIoOgEKz2WyaOHGi4/H8+fNls9l04MAByzJdrl+/foqKisozuzx3cVm/fr1sNpvWr1/vmN1yyy1q1KhRsb+3JB04cEA2m03z5893yfsB7oyiA1jkxx9/VM+ePVWzZk0FBQWpWrVq6tixo15//XWro+FPFixYoGnTplkdwyl3zga4Cxv3ugJcb+PGjbr11ltVo0YNJSQkqHLlyjp8+LC+/fZb7d27V3v27LE64hXZbDZNmDDBsXXEbrfrwoULCgwMlM1mszbc/+vXr5/Wr1+fZyvT+fPn5efnJz8/vwK/zh133KGdO3cWamtVbm6usrOzFRAQIB+fi/+evOWWW5SamqqdO3cW+HWuNZthGMrKypK/v798fX2L7P2Akqjg/7cDKDLPP/+8QkNDtWXLFpUrVy7Pc8eOHbMm1HXw9fUtEb9Qg4KCivX1z58/7yg3xf1eV2Kz2Sx9f8CdsOsKsMDevXsVGxtrKjmSFB4enufxvHnz9Je//EXh4eEKDAxUTEyMZsyYYVovKipKd9xxh9avX6+WLVsqODhYjRs3dhwnsmTJEjVu3FhBQUFq0aKFtm/fnmf9fv36qUyZMtq3b586d+6s0qVLq2rVqnrmmWd0tQ2/zo7RuZTn66+/VuvWrRUUFKRatWrp3XffNa3/ww8/qEOHDgoODlb16tX13HPPad68eQU+7mfp0qVq1KiRgoKC1KhRI3388cdOl7v8GJ2MjAyNGDFCUVFRCgwMVHh4uDp27Kht27ZJurgVZtmyZTp48KBsNptsNpvjuJ9Lx+EsXLhQY8eOVbVq1VSqVCmlp6c7PUbnkq1bt6pt27YKDg5WdHS0Zs6cedXP8s/vd+k1r5Qtv2N0/vOf/6h9+/YqXbq0ypUrpzvvvFO7du3Ks8zEiRNls9m0Z88e9evXT+XKlVNoaKj69++vs2fP5v+HALgptugAFqhZs6Y2bdqknTt3XvUA1RkzZig2NlZ//etf5efnp88++0yPPvqocnNzNWTIkDzL7tmzR/fee68GDx6s+++/Xy+//LK6d++umTNnavTo0Xr00UclSUlJSerVq5d2797t2LUiXdwFdfvtt+vGG2/Uiy++qBUrVmjChAnKycnRM888U+jvc8+ePerZs6cGDhyohIQEzZ07V/369VOLFi0UGxsrSTpy5IhuvfVW2Ww2jRo1SqVLl9Zbb72lwMDAAr3HqlWrdM899ygmJkZJSUk6ceKE+vfvr+rVq1913YcffliLFy/W0KFDFRMToxMnTujrr7/Wrl271Lx5c40ZM0ZpaWn67bff9Morr0iSypQpk+c1nn32WQUEBOiJJ55QVlaWAgIC8n2/U6dOqWvXrurVq5f69OmjDz74QI888ogCAgI0YMCAAn2/lxQk25+tWbNGXbp0Ua1atTRx4kSdO3dOr7/+utq1a6dt27aZDtzu1auXoqOjlZSUpG3btumtt95SeHi4Jk+eXKicgOUMAC63atUqw9fX1/D19TXi4uKMJ5980li5cqWRnZ1tWvbs2bOmWefOnY1atWrlmdWsWdOQZGzcuNExW7lypSHJCA4ONg4ePOiYz5o1y5BkrFu3zjFLSEgwJBnDhg1zzHJzc41u3boZAQEBxvHjxx1zScaECRMcj+fNm2dIMvbv32/Ks2HDBsfs2LFjRmBgoPH44487ZsOGDTNsNpuxfft2x+zEiRNGhQoVTK/pTLNmzYwqVaoYp0+fdsxWrVplSDJq1qyZZ9nLc4eGhhpDhgy54ut369bN9DqGYRjr1q0zJBm1atUy/Rldeu7Pn2+HDh0MScaUKVMcs6ysLKNZs2ZGeHi448/e2WeZ32vml23//v2GJGPevHmO2aX3OXHihGP2/fffGz4+Pkbfvn0dswkTJhiSjAEDBuR5zbvuusuoWLGi6b0Ad8euK8ACHTt21KZNm/TXv/5V33//vV588UV17txZ1apV06effppn2eDgYMd/p6WlKTU1VR06dNC+ffuUlpaWZ9mYmBjFxcU5Hrdp00aS9Je//EU1atQwzfft22fKNnToUMd/22w2DR06VNnZ2VqzZk2hv8+YmBi1b9/e8TgsLEz169fP874rVqxQXFycmjVr5phVqFBB991331Vf//fff9eOHTuUkJCg0NBQx7xjx46KiYm56vrlypXTf//7Xx09erSA35FZQkJCnj+jK/Hz89PgwYMdjwMCAjR48GAdO3ZMW7duveYMV3Ppc+rXr58qVKjgmDdp0kQdO3bUF198YVrn4YcfzvO4ffv2OnHihNLT04stJ1AcKDqARVq1aqUlS5bo1KlT2rx5s0aNGqWMjAz17NlTP/30k2O5b775RvHx8Y7jKsLCwjR69GhJMhWdP5cZSY5f/pGRkU7np06dyjP38fFRrVq18szq1asnSdd0jZzL80hS+fLl87zvwYMHVadOHdNyzmaXO3jwoCSpbt26pufq169/1fVffPFF7dy5U5GRkWrdurUmTpzotPxdSXR0dIGXrVq1qkqXLp1ndj2fb0Fd+pycfSYNGzZUamqqMjMz88wv/7MrX768JPPPDODuKDqAxQICAtSqVStNmjRJM2bM0IULF/Thhx9KunjQ8m233abU1FRNnTpVy5Yt0+rVq/XYY49Junga85/ld+ZTfnOjmK8uYdX7FlSvXr20b98+vf7666patapeeuklxcbGavny5QV+jYJuzSmo/E7Pt9vtRfo+V+Puf3ZAQVF0ADfSsmVLSRd3NUjSZ599pqysLH366acaPHiwunbtqvj4+CL/5XpJbm6uaYvGL7/8Ikmmg1WLSs2aNZ1eN6gg1xKqWbOmJOnXX381Pbd79+4CvX+VKlX06KOPaunSpdq/f78qVqyo559/3vF8UV4X6OjRo6YtJ5d/vpe2nJw+fTrPcpe2yvxZQbNd+pycfSY///yzKlWqZNrSBHgKig5ggXXr1jn9l/GlYyUu7WK49K/qPy+blpamefPmFVu2N954w/HfhmHojTfekL+/v2677bZieb/OnTtr06ZN2rFjh2N28uRJvf/++1ddt0qVKmrWrJneeeedPLvxVq9enWf3nzN2u9206y88PFxVq1ZVVlaWY1a6dGnTctcqJydHs2bNcjzOzs7WrFmzFBYWphYtWkiSateuLUnasGFDnqyzZ882vV5Bs/35c/pzgdq5c6dWrVqlrl27Xuu3BLg9Ti8HLDBs2DCdPXtWd911lxo0aKDs7Gxt3LhRixYtUlRUlPr37y9J6tSpkwICAtS9e3cNHjxYZ86c0Zw5cxQeHu7Y6lOUgoKCtGLFCiUkJKhNmzZavny5li1bptGjRyssLKzI30+SnnzySb333nvq2LGjhg0b5ji9vEaNGjp58uRVt1okJSWpW7duuummmzRgwACdPHlSr7/+umJjY3XmzJl818vIyFD16tXVs2dPNW3aVGXKlNGaNWu0ZcsWTZkyxbFcixYttGjRIiUmJqpVq1YqU6aMunfvfk3fa9WqVTV58mQdOHBA9erV06JFi7Rjxw7Nnj1b/v7+kqTY2FjdeOONGjVqlE6ePKkKFSpo4cKFysnJMb1eYbK99NJL6tKli+Li4jRw4EDH6eWhoaEuuf8XYBkrT/kCvNXy5cuNAQMGGA0aNDDKlCljBAQEGHXq1DGGDRtmpKSk5Fn2008/NZo0aWIEBQUZUVFRxuTJk425c+c6PZ27W7dupveSZDqF+tLpxy+99JJjlpCQYJQuXdrYu3ev0alTJ6NUqVJGRESEMWHCBMNut5tesyCnlzvL06FDB6NDhw55Ztu3bzfat29vBAYGGtWrVzeSkpKM1157zZBkJCcn5/cxOnz00UdGw4YNjcDAQCMmJsZYsmSJkZCQcMXTy7Oysox//vOfRtOmTY2QkBCjdOnSRtOmTY0333wzzzpnzpwx7r33XqNcuXJ5Tlm/dLr3hx9+aMqT3+nlsbGxxnfffWfExcUZQUFBRs2aNY033njDtP7evXuN+Ph4IzAw0IiIiDBGjx5trF692vSa+WVzdnq5YRjGmjVrjHbt2hnBwcFG2bJlje7duxs//fRTnmUunV7+58sJGEb+p70D7o57XQGQdPHKyIsXL77iVhBXGjFihGbNmqUzZ86UiNtLAHBPHKMDwHLnzp3L8/jEiRP617/+pZtuuomSA+C6cIwOAMvFxcXplltuUcOGDZWSkqK3335b6enpGjdunNXRAJRwFB0AluvatasWL16s2bNny2azqXnz5nr77bd18803Wx0NQAnHMToAAMBjcYwOAADwWBQdAADgsbzuGJ3c3FwdPXpUISEhRXppdwAAUHwMw1BGRoaqVq0qH5+Cb6fxuqJz9OhR052cAQBAyXD48GFVr169wMt7XdEJCQmRdPGDKlu2rMVpAABAQaSnpysyMtLxe7ygvK7oXNpdVbZsWYoOAAAlTGEPO+FgZAAA4LEoOgAAwGNRdAAAgMei6AAAAI9F0QEAAB6LogMAADwWRQcAAHgsig4AAPBYFB0AAOCxKDoAAMBjWVp0NmzYoO7du6tq1aqy2WxaunTpVddZv369mjdvrsDAQNWpU0fz588v9pwAAKBksrToZGZmqmnTppo+fXqBlt+/f7+6deumW2+9VTt27NCIESP04IMPauXKlcWc1H0YhiG73Z7nyzAMq2PBg/AzBuBaZJ7L0uHkU46v5NR0qyNJsvimnl26dFGXLl0KvPzMmTMVHR2tKVOmSJIaNmyor7/+Wq+88oo6d+5cXDHdSm5urvbv359nFh0dLV9fX4sSwdPwM4brZRiGjp86ows59gKvs/+3VB06elJ+fp5/RMX/9hzVzl+PKrJKBaujFKnDv5/M87hGlQp6ZWQvi9L8oUTdvXzTpk2Kj4/PM+vcubNGjBiR7zpZWVnKyspyPE5Pd4+GCQBFKSPzvDIyz1/z+r+lnNbeQ8fk71+wXwvLN+xU9oUcVSpfJs887cw5pWWcu+Yc3uTyYoDiUaKKTnJysiIiIvLMIiIilJ6ernPnzik4ONi0TlJSkp5++mlXRQSAa3L2XLZ+SzmlHT8flq+veavGvsOp+vb7fYqqVinP3DAMHTx6wlUxTQ7xyxr5KFsmyOoIkkpY0bkWo0aNUmJiouNxenq6IiMjLUwEwBNlX8jRrwePadtPh7R07Q5FVi4vP7+C7e7b/1tqgd/nwJGCLwtYxSapS/tGVseQVMKKTuXKlZWSkpJnlpKSorJlyzrdmiNJgYGBCgwMdEU8AB7KMAz9lnJaX2/bo8Urt6pGlQry/1OJ2Xv4uGmdw8mnXBnRozSPqaGAApbEkirXMPR7arru6NBYQYH+VscpUjabTXVqhCmiYlmro0gqYUUnLi5OX3zxRZ7Z6tWrFRcXZ1EiAJ4gI/O8Zn3wlZJT0+Trk3e3UXJqms6czcoz8/TdNY3rVVNQwNV/+Waey5KPj03xNzaUzWbL85yvr48a1q6ssqULvvvCZrOZXge4XpYWnTNnzmjPnj2Ox/v379eOHTtUoUIF1ahRQ6NGjdKRI0f07rvvSpIefvhhvfHGG3ryySc1YMAA/ec//9EHH3ygZcuWWfUtAHADGZnnlX0hp1DrZJ7L1tsffa2dvx4tplTXJ7xCiGpcdlaOYUhHjp3S7Tc1cnr8Q40qFRRZufw1v6fNZnN6fBBQklladL777jvdeuutjseXjqVJSEjQ/Pnz9fvvv+vQoUOO56Ojo7Vs2TI99thjevXVV1W9enW99dZbXnNqOYCLvtm+V6/+a62qVArVbynuuYuodHCg2t5QSw2iKxdoecOQqoSFqnZkmPz9PXu3DeBKNsPLrgSWnp6u0NBQpaWlqWxZ99h/WBh2u51rnKBYudPPWEbmeWVl5+jZGRe32pYpHaif9yW7PIczHds2VP2ovCWmXNlSiqldWYEF2O0DoHCu9fd3iTpGB4DnOJmWqe92HpQ9N1eS9O33+5Scmq7wCiFKPXVGx05mWJIrwN9PA+9pp8sPFfH381VM7aqqWK40x5EAJQhFB0CxO3M2S3Z7rn789Yjmf7xRp9LP5rts6qkzLkz2h9EPdVHT+tULfEo4gJKBogOgyGVkntfyr3Zq0fLviv297vxLU1UJC1W9qAjTVXoLolRQAFtoAA9G0QFw3ez2XJ3PvqALOXYNHPtusb7XQ39rL0kqUzpIrRrVVEABb1kAwDvxNwSAa5aVfUH3/vPt636doEB/1apeSeeyLuh0+lndFX+DfH18ZLNJDWtXUbmQYJUODuTUZwCFRtEBUGBZ2ReUciJDu/cna8GyLUo/c+03b6wfXVm9u7RUk3rV2HUEoNhQdABc1b8+/Vafrf/xmtfv3C5WnW+KVdkyQSoTHMh1YgC4DEUHgFMn0s5q/5GTmv/pVpUqVeqatrr8++UHOYYGgKX4GwiApIu7pXb+elT7fjuutz/88ppe49XRvRUU4K/yZUtxPA0At0DRAbzc+awLGp60yHH9msJcLN3HZlPtGmF6tM8tpvsyAYA7oOgAXmz5Vzv11uKvC7VOYr+OalSnqkJDgospFQAUHYoO4CVycuzafSBFR4+d1sxFGwq1bplSgcqxG3pzfB+VCylVTAkBoOhRdAAvYBiGej8+p1DrVK4UopubR+vm5tHcOBZAiUXRATzc4InvFer+UTMm3KuMU8eLMREAuA5FB/BQx05m6JGn3y/w8hMevUNN6leX3W6n6ADwGBQdwMNkX8jRU1OW6NDvJ6+6bKtGURpwTzuFlS/D1YkBeCSKDuBBdvx8WM/OWHbV5T585SH5+HCdGwCej6IDeIjFq7bp38s2X3GZ+ZP6KaR0kIsSAYD1KDpACbf1fwc1afbyKy7TIqamRg/u4qJEAOA+KDpACZSbm6t9h1P11NQlV132vckDFBwU4IJUAOB+KDpACZN9IUd9nnjrqst1v6WJEnrEcZAxAK9G0QFKkIKWnA6t6qnfXW1dkAgA3BtFBygh3liwTuv+u/uKy3Tr0Fh3/qWpKpYr46JUAODeKDpACXDP8JlXfD7A309zn+vLsTgAcBmKDuDmtv7v4BWfXzxtMMfhAEA+uGIY4MYMw7jiqeOUHAC4MrboAG5s/sebnM5fHd1b1SPKuzgNAJQ8bNEB3NT6zbv1+Zc/mOaP/L0DJQcACogtOoAbmrHwS63ZtMvpc7fd2MDFaQCg5GKLDuBm1m/enW/JGTO4K8fkAEAhsEUHcANpGef0wlsr9MuBlHyXqRcVoeYxNVyYCgBKPooOYLGcHLsGjH3niss8lhCvm5rXcVEiAPAc7LoCLLb0P99f8XlfXx9KDgBcI7boABb797LN+T73YM+b1KV9IxemAQDPQtEBLLTgc+cl542xfVQlLNTFaQDA87DrCrDIByu+00ert5nmt93YgJIDAEWEogNYIDk1XYuWf+f0uUE927s4DQB4LooOYIEhzy5wOh/9UBf5+/u6OA0AeC6O0QFc6Mix01q2/kenzz3ap4NaxNZ0cSIA8GwUHcBFfvzliCZO/8zpc13aN9JtNzZ0cSIA8HzsugJcICfHnm/JkS6eRg4AKHps0QGKkWEY+nDl1nwPPJak7rc0cWEiAPAuFB2gGD3z5jL98Mtv+T4/qGd7db4pxoWJAMC7UHSAYpJ5LuuKJeejVx92YRoA8E4cowMUk74j5+X73OJpg12YBAC8F1t0gGLw4y9HnM4nJ96tOjXDXZwGALwXRQcoYt/v/k3PvPm50+coOQDgWhQdoAjNW7JRn3/5g9Pnpo3q7eI0AACO0QGKiGEY+ZacPt1aK7JyeRcnAgBQdIAi8t3/Djqd14+urJ6dmrs4DQBAYtcVUGRmLtxgmvXu0lK9bm9pQRoAgMQWHaDInM44a5pRcgDAWhQdoAjcM3ym1REAAE5QdIDr9L89R53O3xjbx8VJAACXo+gA12n86586nVcJC3VxEgDA5Sg6wHXIvpDjdP7Ws31dnAQA4IzlRWf69OmKiopSUFCQ2rRpo82bN19x+WnTpql+/foKDg5WZGSkHnvsMZ0/f95FaYE/5OTYNfS5f5vmDWpVVvmypSxIBAC4nKVFZ9GiRUpMTNSECRO0bds2NW3aVJ07d9axY8ecLr9gwQKNHDlSEyZM0K5du/T2229r0aJFGj16tIuTw9tlX8hR78fn6MTpTNNzox/qYkEiAIAzlhadqVOnatCgQerfv79iYmI0c+ZMlSpVSnPnznW6/MaNG9WuXTvde++9ioqKUqdOndSnT5+rbgUCilJubq76PPFWvs+XDg50YRoAwJVYVnSys7O1detWxcfH/xHGx0fx8fHatGmT03Xatm2rrVu3OorNvn379MUXX6hr1675vk9WVpbS09PzfAHXY9D49/J9LrFfRxcmAQBcjWVXRk5NTZXdbldERESeeUREhH7++Wen69x7771KTU3VTTfdJMMwlJOTo4cffviKu66SkpL09NNPF2l2eK+cHLvTCwNK0l23NVO7G2q7OBEA4EosPxi5MNavX69JkybpzTff1LZt27RkyRItW7ZMzz77bL7rjBo1SmlpaY6vw4cPuzAxPE3vx+c4nS98eZDu/+uNLk4DALgay7boVKpUSb6+vkpJSckzT0lJUeXKlZ2uM27cOD3wwAN68MEHJUmNGzdWZmamHnroIY0ZM0Y+PubeFhgYqMBAjpnA9ctvS86MCffJ39/XxWkAAAVh2RadgIAAtWjRQmvXrnXMcnNztXbtWsXFxTld5+zZs6Yy4+t78ReMYRjFFxZeLyfHroFj33X6XHiFEBenAQAUlKV3L09MTFRCQoJatmyp1q1ba9q0acrMzFT//v0lSX379lW1atWUlJQkSerevbumTp2qG264QW3atNGePXs0btw4de/e3VF4gOKQ3y6rD6Y+5OIkAIDCsLTo9O7dW8ePH9f48eOVnJysZs2aacWKFY4DlA8dOpRnC87YsWNls9k0duxYHTlyRGFhYerevbuef/55q74FeIGUE87P1IuqVkm+viXqMDcA8Do2w8v2+aSnpys0NFRpaWkqW7as1XEKzW63a//+/Xlm0dHRbNEqRvndmfyjVx92cRLX4GcMgDu61t/f/HMUuAaeWnIAwNNQdIAreOAp81W6y5TiLD4AKCkoOkA+MjLP6+z5bNP8ueE9XB8GAHBNKDpAPvqNnu90Hlm5vGuDAACuGUUHcGLukm+czt+bPMDFSQAA14OiA1zmwgW7ln35o2keVa2SgoMCLEgEALhWFB3gMouWb3E6f3bYX12cBABwvSg6wGU+XrvDNJs2qrdKBbM1BwBKGooO8Cfbdzm/uz0HIANAyUTRAf7kuZnLTLPH+3e0IAkAoChQdID/l9/dUNo2q+3iJACAokLRAf7fgSMnTLMm9apbkAQAUFQoOsD/e+KlxabZmMFdLEgCACgqFB1A0onTZ5zO/fy4YzcAlGQUHUDSQxPeM81qVKlgQRIAQFGi6MDr5ebmOp2/kHiXi5MAAIoaRQde72+PzTbNwsqHKDDA34I0AICiRNGBV7tn+Eyn89fG9HZxEgBAcaDowGt9tHpbvs8F+Pu5MAkAoLhQdOCV0jLOacHnm50+N+/5BBenAQAUF4oOvNKAse84nb85/l6VLRPs4jQAgOJC0YHXye/GnUmP3aWIimVdnAYAUJwoOvAq2RdynN64MzDAX/WiIixIBAAoThQdeJXX3lvndL7gpYEuTgIAcAWKDrzKph17TbOEHnEWJAEAuAJFB14jK/uC0/lfb23q4iQAAFeh6MBrfLRqu2k2ZnBXC5IAAFyFogOv4ewCgc1jaliQBADgKhQdeIWf9v5udQQAgAUoOvB4J9MyNe61T0zzuGa1LUgDAHAlig483lNTljidP9b3NhcnAQC4GkUHHu9kWqZp1rldrHx9+fEHAE/H3/TwaIZhmGY2SQ/1au/6MAAAl6PowGMZhqGeI2aZ5q+N+bsFaQAAVqDowGM5KzmSVCUs1MVJAABWoejAI+05eCzf52w2mwuTAACsRNGBxzn0+0k9NdX5mVbvvtDfxWkAAFbyszoAUNQee+EDp/N/v/ygAvz5kQcAb8IWHXiUzHNZTue9u7Sk5ACAF6LowKP0HTnPNAsNCVav21takAYAYDWKDjze28/2tToCAMAiFB14jAsX7E7nnGUFAN6LogOP8fcn5phmTw/tbkESAIC7oOjAIySnpjudx9Su4uIkAAB3QtGBRxjy7ALTzCbJx4cfcQDwZvwWQIlnt+c6nS+a+pCLkwAA3A1FByXekjXbTbMOrerJ15cfbwDwdvwmQIm38Istptmw+261IAkAwN1QdFCivffpt07nnFIOAJAoOijBjp3M0Mdrd5jmndvFuj4MAMAtUXRQYg2ftMjpfOA97VycBADgrig6KLGyL+SYZlOf+hsHIQMAHPiNAI9Ss2pFqyMAANwIRQcl0r+dnGnVs3MLC5IAANwZRQcl0uKVW02zezreYEESAIA7o+jAYwT4+1kdAQDgZiwvOtOnT1dUVJSCgoLUpk0bbd68+YrLnz59WkOGDFGVKlUUGBioevXq6YsvvnBRWriDg0dPmmb33dHGgiQAAHdn6T+BFy1apMTERM2cOVNt2rTRtGnT1LlzZ+3evVvh4eGm5bOzs9WxY0eFh4dr8eLFqlatmg4ePKhy5cq5Pjwss+LrnabZ7Tdx7RwAgJmlRWfq1KkaNGiQ+vfvL0maOXOmli1bprlz52rkyJGm5efOnauTJ09q48aN8vf3lyRFRUW5MjLcwKpvfjLNSgUHWJAEAODuLNt1lZ2dra1btyo+Pv6PMD4+io+P16ZNm5yu8+mnnyouLk5DhgxRRESEGjVqpEmTJslut+f7PllZWUpPT8/zhZJr0459VkcAAJQglhWd1NRU2e12RURE5JlHREQoOTnZ6Tr79u3T4sWLZbfb9cUXX2jcuHGaMmWKnnvuuXzfJykpSaGhoY6vyMjIIv0+4Fovz1tlmv311qYWJAEAlASWH4xcGLm5uQoPD9fs2bPVokUL9e7dW2PGjNHMmTPzXWfUqFFKS0tzfB0+fNiFieEKfe+80eoIAAA3ZdkxOpUqVZKvr69SUlLyzFNSUlS5cmWn61SpUkX+/v7y9fV1zBo2bKjk5GRlZ2crIMB8nEZgYKACAwOLNjwssfnHA6bZDQ0juVM5ACBflm3RCQgIUIsWLbR27VrHLDc3V2vXrlVcXJzTddq1a6c9e/YoNzfXMfvll19UpUoVpyUHnuN81gVNfmuFaZ7Qo60FaQAAJYWlu64SExM1Z84cvfPOO9q1a5ceeeQRZWZmOs7C6tu3r0aNGuVY/pFHHtHJkyc1fPhw/fLLL1q2bJkmTZqkIUOGWPUtwEUmzV7udB5ZubyLkwAAShJLTy/v3bu3jh8/rvHjxys5OVnNmjXTihUrHAcoHzp0SD4+f3SxyMhIrVy5Uo899piaNGmiatWqafjw4Xrqqaes+hbgIv/bc9Q0G9SzvQVJAAAlieXXzB86dKiGDh3q9Ln169ebZnFxcfr222+LORXcycYde53Ob2/PRQIBAFdWos66gneaMm+1afbi4/dYkAQAUNJQdFAi1a4RZnUEAEAJQNGBWzufdcE0a9e8jgVJAAAlEUUHbu3DlVtNs349nF9+AACAy1F04Nac3cCzQmhpC5IAAEoiig7c2tnz2VZHAACUYBQduK3Mc1mmWc2qFS1IAgAoqSg6cFt9R84zze69o7UFSQAAJRVFB27pcPIpp/MWMTVcnAQAUJJRdOCW3v1kk2lWs2pF7lQOACgUig7c0rafDplmU57saUESAEBJRtGB2zEMw+mcrTkAgMKi6MDtPDVliWnW/662FiQBAJR0FB24lYNHT2jv4eOmebcOjS1IAwAo6Sg6cCuJkz90Ome3FQDgWhSq6PTt21cZGRmOx99//70uXDDfdBG4FufyuQry1Kd6uTgJAMBTFKrovP/++zp37pzjcfv27XX48OEiDwXvNHH656bZzS3rqmbVChakAQB4gkIVncvPhsnv7BjgWuw5dMw0+8f9f7EgCQDAU3CMDtwCp5QDAIqDX2FX+Omnn5ScnCzp4i+nn3/+WWfOnMmzTJMmTYomHbzGvf982zTr0437WgEArk+hi85tt92W51/fd9xxh6SL//I2DEM2m012u73oEsLjGYah7As5pnmPvzS1IA0AwJMUqujs37+/uHLAiy38YovTuZ+fr4uTAAA8TaGKTs2aNYsrB7zY4lXbTLP5k/q5PggAwOMUeteVJP3666/65JNPdODAAdlsNkVHR6tHjx6qVatWUeeDlwopHWR1BACAByh00UlKStL48eOVm5ur8PBwGYah48ePa+TIkZo0aZKeeOKJ4sgJD/WvT781zQbc3c6CJAAAT1So08vXrVunsWPHasyYMUpNTdXvv/+u5ORkR9EZOXKkNmzYUFxZ4WFycuxaunaHaX5Dw0jXhwEAeKRCbdGZOXOmHnzwQU2cODHPvEKFCnrmmWeUnJysGTNm6Oabby7KjPBQvR+f43ReNbyca4MAADxWobbobN68WQ888EC+zz/wwAP69lvzrgjgcvldIHDqU39zcRIAgCcrVNFJSUlRVFRUvs9HR0c7LiYIXMmufc5/TmpWrejiJAAAT1aoonP+/HkFBATk+7y/v7+ys53fgRr4s6ffNN/A8/0XB1qQBADgyQp91tVbb72lMmXKOH0uIyPjugPBO+TkmK+eHRTob0ESAIAnK1TRqVGjhubMcX4A6Z+XAa6Eu94DAFylUEXnwIEDxRQD3mTA2HdNs753xlmQBADg6Qp1jM5//vMfxcTEKD093fRcWlqaYmNj9dVXXxVZOHgewzCUfuacaR7XjKtqAwCKXqGKzrRp0zRo0CCVLVvW9FxoaKgGDx6sqVOnFlk4eJ5fDqQ4nYdXCHFxEgCANyhU0fn+++91++235/t8p06dtHXr1usOBc+19X+HTLP3Jg+wIAkAwBsU+jo6/v75nxnj5+en48ePX3coeK6PVpvvVB4clP8lCwAAuB6FKjrVqlXTzp07833+hx9+UJUqVa47FLyHv5+v1REAAB6sUEWna9euGjdunM6fP2967ty5c5owYYLuuOOOIgsHz3L02GnT7JbW9VwfBADgNQp1evnYsWO1ZMkS1atXT0OHDlX9+vUlST///LOmT58uu92uMWPGFEtQlHzDnl9omvXr0daCJAAAb1GoohMREaGNGzfqkUce0ahRoxwXfrPZbOrcubOmT5+uiIiIYgkKz8TVkAEAxanQt4CoWbOmvvjiC506dUp79uyRYRiqW7euypcvXxz54CE+W/eDaUbJAQAUt0IXnUvKly+vVq1aFWUWeLD5SzeaZrMm3m9BEgCANynUwchAUSpTKtDqCAAAD0fRQbE74uRsq7Y31HZ9EACA16HooNj9w8nZVn26stsTAFD8KDqwRNXwclZHAAB4AYoOitXeQ+ZbgjSpV92CJAAAb0TRQbF6cspHptmYwV0sSAIA8EYUHRSbzHNZTud+3N8KAOAiFB0Um3++ZN6a07AWN30FALgORQfFJuVEumk2/tFuFiQBAHgrig6Kxbnz2U7nAf7XfDFuAAAKjaKDYvHQhPdMs9fG/N2CJAAAb0bRQbE462SLTtWwUAuSAAC8GUUHRc5uz3U6t9lsLk4CAPB2blF0pk+frqioKAUFBalNmzbavHlzgdZbuHChbDabevToUbwBUShHj6eZZg/2vMmCJAAAb2d50Vm0aJESExM1YcIEbdu2TU2bNlXnzp117NixK6534MABPfHEE2rfvr2LkqKgFi4zF9Xbb4q1IAkAwNtZXnSmTp2qQYMGqX///oqJidHMmTNVqlQpzZ07N9917Ha77rvvPj399NOqVauWC9OiIL79Yb9pxm4rAIAVLC062dnZ2rp1q+Lj4x0zHx8fxcfHa9OmTfmu98wzzyg8PFwDBw50RUwUQkbmeasjAADgYOlFTVJTU2W32xUREZFnHhERoZ9//tnpOl9//bXefvtt7dixo0DvkZWVpaysP25FkJ5uvogdis7cJd+YZk89eLsFSQAAcINdV4WRkZGhBx54QHPmzFGlSpUKtE5SUpJCQ0MdX5GRkcWc0rv9sPuIada6cZTrgwAAIIu36FSqVEm+vr5KSUnJM09JSVHlypVNy+/du1cHDhxQ9+7dHbPc3IunMvv5+Wn37t2qXbt2nnVGjRqlxMREx+P09HTKTjE6nXHW6ggAADhYWnQCAgLUokULrV271nGKeG5urtauXauhQ4ealm/QoIF+/PHHPLOxY8cqIyNDr776qtMCExgYqMDAwGLJDzObJONPjxvUMhdWAABcxfIbDyUmJiohIUEtW7ZU69atNW3aNGVmZqp///6SpL59+6patWpKSkpSUFCQGjVqlGf9cuXKSZJpDmsYlz1uHlPDkhwAAEhuUHR69+6t48ePa/z48UpOTlazZs20YsUKxwHKhw4dko9PiTqUyGs5u1t5nRrhFiQBAOAim2EYl/8j3KOlp6crNDRUaWlpKlu2rNVxCs1ut2v//rzXqYmOjpavr69Fif4weOJ7Sj11Js/s3y8/yB3LSxh3/hkD4L2u9fc3m0pQJLIv5JhKjiRKDgDAUhQdFImX5q4yzcqWCbYgCQAAf6DooEhs++mQaTblyZ4WJAEA4A8UHVy381kXnM4rhJZ2cRIAAPKi6OC6vfvJt6bZ3fE3WJAEAIC8KDq4biu/+Z9p1u2WxhYkAQAgL4oOikW5kFJWRwAAgKKD63Phgt0069wu1oIkAACYUXRwXX45mGKade3A7TgAAO6BooPr8sKcFaZZtfByrg8CAIATFB1cl7Pns00zm81mQRIAAMwoOrhmJ9MyTbPQEK6GDABwHxQdXLNB4/9lmo0adLsFSQAAcI6igyJVt2aE1REAAHCg6KDI+Pv5Wh0BAIA8KDq4JgePnjDNHkuItyAJAAD5o+jgmiz78kfT7IaGkRYkAQAgfxQdXJNde383zQL8/SxIAgBA/ig6uCZHj6dZHQEAgKui6KBI3NqmvtURAAAwoeig0NZ+u8s0qxMZbkESAACujKKDQnvz31+aZre0rmdBEgAAroyigyIRFOhvdQQAAEwoOiiU3fuTTbNuHRpbkAQAgKuj6KBQ5iz+2jS7v3sbC5IAAHB1FB0Uyv7fUk0zrp8DAHBXFB0UWE6O3eoIAAAUCkUHBbbPydacx/t3tCAJAAAFQ9FBgY165WPTLK5pLQuSAABQMBQdFMhHq7c5ndtsNhcnAQCg4Cg6KJAFn282zbhbOQDA3VF0cM3GPtzN6ggAAFwRRQdX9c32vaZZ15sbWZAEAIDCoejgqqbOX22a9e7SyoIkAAAUDkUHV2QYhtN5mVKBLk4CAEDhUXRwRbM+2GCa3XVbM9cHAQDgGlB0cEWrN+4yzXp2bmFBEgAACo+ig3ydTMt0Og8K9HdxEgAArg1FB/l6/b11ptmjfTpYkAQAgGtD0UG+fvjlN9PsthsbWpAEAIBrQ9EBAAAei6IDp77fbd6acyM38AQAlDAUHTj1zJufm2btW9SxIAkAANeOogOTd5Zucjpv0yTaxUkAALg+FB2YfLrue9Psr7c2lc1msyANAADXjqKDPHJy7E7nCT3iXJwEAIDrR9FBHj/+etQ0e/HxeyxIAgDA9aPoII/3PvuvaVYrspIFSQAAuH4UHeRx4EiqacaxOQCAkoqigyuKqV3F6ggAAFwzig4cvv1+n2l2T6fmFiQBAKBoUHTg8NLcVaZZ47rVLEgCAEDRoOjginx9+REBAJRc/BaDJCnzXJZp1iKmpgVJAAAoOhQdSJJWfP0/06xvjxstSAIAQNGh6ECStODzzaZZ9YjyFiQBAKDouEXRmT59uqKiohQUFKQ2bdpo82bzL91L5syZo/bt26t8+fIqX7684uPjr7g8AADwXpYXnUWLFikxMVETJkzQtm3b1LRpU3Xu3FnHjh1zuvz69evVp08frVu3Tps2bVJkZKQ6deqkI0eOuDi559h3+LhpVi8qwoIkAAAULcuLztSpUzVo0CD1799fMTExmjlzpkqVKqW5c+c6Xf7999/Xo48+qmbNmqlBgwZ66623lJubq7Vr17o4uef4748HTLMxg7u6PggAAEXM0qKTnZ2trVu3Kj4+3jHz8fFRfHy8Nm3aVKDXOHv2rC5cuKAKFSoUV0yPt/2nQ6ZZmVKBFiQBAKBo+Vn55qmpqbLb7YqIyLubJCIiQj///HOBXuOpp55S1apV85SlP8vKylJW1h+nTqenp197YA+118muKwAAPIHlu66uxwsvvKCFCxfq448/VlBQkNNlkpKSFBoa6viKjIx0ccqSp0m96lZHAACgSFhadCpVqiRfX1+lpKTkmaekpKhy5cpXXPfll1/WCy+8oFWrVqlJkyb5Ljdq1CilpaU5vg4fPlwk2T1FRuZ506xZQ8ogAMAzWFp0AgIC1KJFizwHEl86sDguLi7f9V588UU9++yzWrFihVq2bHnF9wgMDFTZsmXzfOEP0xesN83at6jj+iAAABQDS4/RkaTExEQlJCSoZcuWat26taZNm6bMzEz1799fktS3b19Vq1ZNSUlJkqTJkydr/PjxWrBggaKiopScnCxJKlOmjMqUKWPZ91FSbdl5wDSrEFra9UEAACgGlhed3r176/jx4xo/frySk5PVrFkzrVixwnGA8qFDh+Tj88eGpxkzZig7O1s9e/bM8zoTJkzQxIkTXRkdAAC4OcuLjiQNHTpUQ4cOdfrc+vXr8zw+cOBA8QfyErv2/m6a9b0z/12GAACUNCX6rCtcn7GvfWKadW3fyIIkAAAUD4qOl8o8l+V07u/v6+IkAAAUH4qOl/rfHvNuq4a1qliQBACA4kPR8VKzP9hgmk0ccocFSQAAKD4UHS908OgJnUo/a5r7+bHbCgDgWSg6XsYwDCVO/tA0rxUZZkEaAACKF0XHy7w8d5XTee8uV77CNAAAJRFFx8t8+8N+06xKWKhaxNSwIA0AAMXLLS4YCGu9MbaP1REAACgWbNHxIiu++p9pdu8drS1IAgCAa1B0vMicxV+ZZnfe2tSCJAAAuAZFx0ucTMt0OueUcgCAJ6PoeIkX315pmnEDTwCAp6PoeAHDMPTrwWOmebebuYEnAMCzUXS8wKfrfnA6Z7cVAMDTUXS8wLufbDLNnnrwdguSAADgWhQdD/dZPltzWjeOcm0QAAAsQNHxYJnnsjR/6UbTvGfnFhakAQDA9Sg6HqzvyHlO5z07NndxEgAArEHR8VDZF3KczhP7dZS/PwchAwC8A0XHQz32wgemWfmypdTuhtoWpAEAwBoUHQ+VnJpumr31bF8LkgAAYB2Kjgd6buYyqyMAAOAWKDoeaPuuw6bZxCHdLUgCAIC1KDoeJmGU8zOtGter5uIkAABYj6LjQXJzc3XmbJZpPmvi/RakAQDAehQdD7J64y6n80rly7g4CQAA7oGi40Fmf/iVafZOUn8LkgAA4B4oOh6uTKlAqyMAAGAZio6HOHj0hGnWrnkdC5IAAOA+KDoeYvzrn5pmA+5ua0ESAADcB0XHAxw9dtrp2VblQkpZkAYAAPdB0fEAo15ZapqVLRPs+iAAALgZio4HOHs+2zR7bvidFiQBAMC9UHRKuN0Hjptm/n6+qhZezvVhAABwMxSdEu71hRtNs2H3/8WCJAAAuB+KTglmGIbTedtmtVycBAAA90TRKcGGTTafUn5L6/qy2WwWpAEAwP1QdEqo5NQMp/NBPW9ycRIAANwXRaeEeu6t/zidBwX6uzgJAADui6JTAr387gan80VTBrk4CQAA7o2iU8IYhqEDR0+Z5tHVK8nPz9eCRAAAuC+KTgnz3mf/dTp/8fG7XZwEAAD3R9EpQez2XH267gfT/PUxveXjwx8lAACX47djCTJ/qfnigJJUuVKoi5MAAFAyUHRKkC827DTN/taxsQVJAAAoGSg6JUBW9gXdM3ym0+c6tOAqyAAA5MfP6gDIn2EY+mLDTs1d8o3T54f9va2LEwEAULJQdNzYwHHvKi3jXL7P148Kc2EaAABKHnZduakFn2++Ysl59tFOLkwDAEDJxBYdN2S35+qj1dvyff6Vx++Qvz8XBwQA4GooOm7EMAzd9+RcZWVfcPr84F436y9t6unAgQOuDQYAQAlF0XEThmGo54hZ+T4//tE71LR+ddntdhemAgCgZKPoWGzh8i36cMXWqy7XtH51F6QBAMCzUHQslN+1cS73rxcGFHMSAAA8E0XHAt9s36up81dfdbn3XxyooEB/FyQCAMAzUXRc5HTGWaWePKN3Ptmkn/b+ftXlF08bLJvN5oJkAAB4LrcoOtOnT9dLL72k5ORkNW3aVK+//rpat26d7/Iffvihxo0bpwMHDqhu3bqaPHmyunbt6sLEV5aWcU5pZ84pOzunwMVGkh7ufbNaN45WaEhwMScEAMA7WF50Fi1apMTERM2cOVNt2rTRtGnT1LlzZ+3evVvh4eGm5Tdu3Kg+ffooKSlJd9xxhxYsWKAePXpo27ZtatSokQXfwUV7Dh7TU1OXXPP6777QX6WDA4swEQAAsBmGYVgZoE2bNmrVqpXeeOMNSVJubq4iIyM1bNgwjRw50rR87969lZmZqc8//9wxu/HGG9WsWTPNnHn1g3vT09MVGhqqtLQ0lS1btki+h5wcu3o/Puea1o2Pa6hH/t6hwMvb7Xbt378/zyw6Olq+vlxAEEWDnzEA7uhaf39beguI7Oxsbd26VfHx8Y6Zj4+P4uPjtWnTJqfrbNq0Kc/yktS5c+d8l8/KylJ6enqer6K277fUa1rvqQdvL1TJAQAAhWPprqvU1FTZ7XZFRETkmUdEROjnn392uk5ycrLT5ZOTk50un5SUpKeffrpoAufjfJbzKxlfbnLi3apT07w7DgAAFA/Lj9EpbqNGjVJiYqLjcXp6uiIjI4v9fV8b83dVCy9X7O8DAADyZ2nRqVSpknx9fZWSkpJnnpKSosqVKztdp3LlyoVaPjAwUIGBxXuQb+N61bR42uA8s+I6NdzHx0fR0dGmGVBU+BkD4Eks/dsrICBALVq00Nq1ax2z3NxcrV27VnFxcU7XiYuLy7O8JK1evTrf5V3BZrOZvorzvXx9ffN8cb0dFCV+xgB4Est3XSUmJiohIUEtW7ZU69atNW3aNGVmZqp///6SpL59+6patWpKSkqSJA0fPlwdOnTQlClT1K1bNy1cuFDfffedZs+ebeW3AQAA3JDlRad37946fvy4xo8fr+TkZDVr1kwrVqxwHHB86NChPJvN27ZtqwULFmjs2LEaPXq06tatq6VLl1p6DR0AAOCeLL+OjqsVx3V0AABA8SqR19EBAAAoThQdAADgsSg6AADAY1F0AACAx6LoAAAAj0XRAQAAHouiAwAAPBZFBwAAeCyKDgAA8FiW3wLC1S5dCDo9Pd3iJAAAoKAu/d4u7A0dvK7oZGRkSJIiIyMtTgIAAAorIyNDoaGhBV7e6+51lZubq6NHjyokJEQ2m83l75+enq7IyEgdPnzYq++1xedwEZ/DRXwOF/E5/IHP4iI+h4sufQ4//fST6tevn+dm31fjdVt0fHx8VL16datjqGzZsl79Q3sJn8NFfA4X8TlcxOfwBz6Li/gcLqpWrVqhSo7EwcgAAMCDUXQAAIDHoui4WGBgoCZMmKDAwECro1iKz+EiPoeL+Bwu4nP4A5/FRXwOF13P5+B1ByMDAADvwRYdAADgsSg6AADAY1F0AACAx6LoAAAAj0XRcZENGzaoe/fuqlq1qmw2m5YuXWp1JJdLSkpSq1atFBISovDwcPXo0UO7d++2OpYlZsyYoSZNmjguAhYXF6fly5dbHctSL7zwgmw2m0aMGGF1FJebOHGibDZbnq8GDRpYHcsSR44c0f3336+KFSsqODhYjRs31nfffWd1LJeKiooy/TzYbDYNGTLE6mguZbfbNW7cOEVHRys4OFi1a9fWs88+y72u3FVmZqaaNm2qAQMG6O6777Y6jiW+/PJLDRkyRK1atVJOTo5Gjx6tTp066aefflLp0qWtjudS1atX1wsvvKC6devKMAy98847uvPOO7V9+3bFxsZaHc/ltmzZolmzZqlJkyZWR7FMbGys1qxZ43js5+d9fz2fOnVK7dq106233qrly5crLCxMv/76q8qXL291NJfasmWL7Ha74/HOnTvVsWNH/e1vf7MwletNnjxZM2bM0DvvvKPY2Fh999136t+/v0JDQ/WPf/yjwK/jff8nWaRLly7q0qWL1TEstWLFijyP58+fr/DwcG3dulU333yzRams0b179zyPn3/+ec2YMUPffvut1xWdM2fO6L777tOcOXP03HPPWR3HMn5+fqpcubLVMSw1efJkRUZGat68eY5ZdHS0hYmsERYWlufxCy+8oNq1a6tDhw4WJbLGxo0bdeedd6pbt26SLm7p+ve//63NmzcX6nXYdQXLpKWlSZIqVKhgcRJr2e12LVy4UJmZmYqLi7M6jssNGTJE3bp1U3x8vNVRLPXrr7+qatWqqlWrlu677z4dOnTI6kgu9+mnn6ply5b629/+pvDwcN1www2aM2eO1bEslZ2drffee08DBgyw5EbUVmrbtq3Wrl2rX375RZL0/fff6+uvvy70RgO26MASubm5GjFihNq1a6dGjRpZHccSP/74o+Li4nT+/HmVKVNGH3/8sWJiYqyO5VILFy7Utm3btGXLFqujWKpNmzaaP3++6tevr99//11PP/202rdvr507dyokJMTqeC6zb98+zZgxQ4mJiRo9erS2bNmif/zjHwoICFBCQoLV8SyxdOlSnT59Wv369bM6isuNHDlS6enpatCggXx9fWW32/X888/rvvvuK9TrUHRgiSFDhmjnzp36+uuvrY5imfr162vHjh1KS0vT4sWLlZCQoC+//NJrys7hw4c1fPhwrV69WkFBQVbHsdSf/4XapEkTtWnTRjVr1tQHH3yggQMHWpjMtXJzc9WyZUtNmjRJknTDDTdo586dmjlzptcWnbfffltdunRR1apVrY7ich988IHef/99LViwQLGxsdqxY4dGjBihqlWrFurngaIDlxs6dKg+//xzbdiwQdWrV7c6jmUCAgJUp04dSVKLFi20ZcsWvfrqq5o1a5bFyVxj69atOnbsmJo3b+6Y2e12bdiwQW+88YaysrLk6+trYULrlCtXTvXq1dOePXusjuJSVapUMRX9hg0b6qOPPrIokbUOHjyoNWvWaMmSJVZHscQ///lPjRw5Un//+98lSY0bN9bBgweVlJRE0YF7MgxDw4YN08cff6z169d75UGGV5Kbm6usrCyrY7jMbbfdph9//DHPrH///mrQoIGeeuopry050sUDtPfu3asHHnjA6igu1a5dO9MlJ3755RfVrFnTokTWmjdvnsLDwx0H43qbs2fPyscn76HEvr6+ys3NLdTrUHRc5MyZM3n+dbZ//37t2LFDFSpUUI0aNSxM5jpDhgzRggUL9MknnygkJETJycmSpNDQUAUHB1uczrVGjRqlLl26qEaNGsrIyNCCBQu0fv16rVy50upoLhMSEmI6Pqt06dKqWLGi1x239cQTT6h79+6qWbOmjh49qgkTJsjX11d9+vSxOppLPfbYY2rbtq0mTZqkXr16afPmzZo9e7Zmz55tdTSXy83N1bx585SQkOCVlxqQLp6d+vzzz6tGjRqKjY3V9u3bNXXqVA0YMKBwL2TAJdatW2dIMn0lJCRYHc1lnH3/kox58+ZZHc3lBgwYYNSsWdMICAgwwsLCjNtuu81YtWqV1bEs16FDB2P48OFWx3C53r17G1WqVDECAgKMatWqGb179zb27NljdSxLfPbZZ0ajRo2MwMBAo0GDBsbs2bOtjmSJlStXGpKM3bt3Wx3FMunp6cbw4cONGjVqGEFBQUatWrWMMWPGGFlZWYV6HZthFPISgwAAACUE19EBAAAei6IDAAA8FkUHAAB4LIoOAADwWBQdAADgsSg6AADAY1F0AACAx6LoACgx+vXrpx49elgdA0AJwgUDAbidAwcOKDo6Wtu3b1ezZs0c87S0NBmGoXLlyhXr+/fr10+nT5/W0qVLi/V9ABQ/77yBBoASKTQ01OoIAEoYdl0BKDa5ublKSkpSdHS0goOD1bRpUy1evFiSdOrUKd13330KCwtTcHCw6tatq3nz5kmS4872N9xwg2w2m2655RZJ5l1Xt9xyi4YNG6YRI0aofPnyioiI0Jw5c5SZman+/fsrJCREderU0fLlyx3r2O12DRw40JGpfv36evXVVx3PT5w4Ue+8844++eQT2Ww22Ww2rV+/XpJ0+PBh9erVS+XKlVOFChV055136sCBA8X3AQK4bmzRAVBskpKS9N5772nmzJmqW7euNmzYoPvvv19hYWH68MMP9dNPP2n58uWqVKmS9uzZo3PnzkmSNm/erNatW2vNmjWKjY1VQEBAvu/xzjvv6Mknn9TmzZu1aNEiPfLII/r444911113afTo0XrllVf0wAMP6NChQypVqpRyc3NVvXp1ffjhh6pYsaI2btyohx56SFWqVFGvXr30xBNPaNeuXUpPT3cUrwoVKujChQvq3Lmz4uLi9NVXX8nPz0/PPfecbr/9dv3www9XzAjAOhyjA6BYZGVlqUKFClqzZo3i4uIc8wcffFBnz57VmTNnVKlSJc2dO9e0bn7H6Fx+7Mwtt9wiu92ur776StLFrTWhoaG6++679e6770qSkpOTVaVKFW3atEk33nij06xDhw5VcnKyY2uTs2N03nvvPT333HPatWuXbDabJCk7O1vlypXT0qVL1alTp2v+rAAUH7boACgWe/bs0dmzZ9WxY8c88+zsbN1www2aOHGi7rnnHm3btk2dOnVSjx491LZt20K/T5MmTRz/7evrq4oVK6px48aOWUREhCTp2LFjjtn06dM1d+5cHTp0SOfOnVN2dnaeQuXM999/rz179igkJCTP/Pz589q7d2+hcwNwDYoOgGJx5swZSdKyZctUrVq1PM8FBgYqMjJSBw8e1BdffKHVq1frtttu05AhQ/Tyyy8X6n38/f3zPLbZbHlml7a+5ObmSpIWLlyoJ554QlOmTFFcXJxCQkL00ksv6b///e9Vv58WLVro/fffNz0XFhZWqMwAXIeiA6BYxMTEKDAwUIcOHVKHDh2cLhMWFqaEhAQlJCSoffv2+uc//6mXX37ZcbyL3W4v8lzffPON2rZtq0cffdQxu3yLTEBAgOm9mzdvrkWLFik8PFxly5Yt8lwAigdFB0CxCAkJ0RNPPKHHHntMubm5uummm5SWlqZvvvlGZcuW1d69e9WiRQvFxsYqKytLn3/+uRo2bChJCg8PV3BwsFasWKHq1asrKCioyE4tr1u3rt59912tXLlS0dHR+te//qUtW7Y4zvSSpKioKK1cuVK7d+9WxYoVFRoaqvvuu08vvfSS7rzzTj3zzDOqXr26Dh48qCVLlujJJ59U9erViyQfgKLF6eUAis2zzz6rcePGKSkpSQ0bNtTtt9+uZcuWKTo6WgEBARo1apSaNGmim2++Wb6+vlq4cKEkyc/PT6+99ppmzZqlqlWr6s477yyyTIMHD9bdd9+t3r17q02bNjpx4kSerTuSNGjQINWvX18tW7ZUWFiYvvnmG5UqVUobNmxQjRo1dPfdd6thw4YaOHCgzp8/zxYewI1x1hUAAPBYbNEBAAAei6IDAAA8FkUHAAB4LIoOAADwWBQdAADgsSg6AADAY1F0AACAx6LoAAAAj0XRAQAAHouiAwAAPBZFBwAAeCyKDgAA8Fj/ByIY3Z4KfBurAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "def SimulateSample(lam=2, n=10, iters=1000):\n",
    "    \"\"\"Sampling distribution of L as an estimator of exponential parameter.\n",
    "\n",
    "    lam: parameter of an exponential distribution\n",
    "    n: sample size\n",
    "    iters: number of iterations\n",
    "    \"\"\"\n",
    "    def VertLine(x, y=1):\n",
    "        thinkplot.Plot([x, x], [0, y], color='0.8', linewidth=3)\n",
    "\n",
    "    estimates = []\n",
    "    for _ in range(iters):\n",
    "        xs = np.random.exponential(1.0/lam, n)\n",
    "        lamhat = 1.0 / np.mean(xs)\n",
    "        estimates.append(lamhat)\n",
    "\n",
    "    stderr = RMSE(estimates, lam)\n",
    "    print('standard error', stderr)\n",
    "\n",
    "    cdf = thinkstats2.Cdf(estimates)\n",
    "    ci = cdf.Percentile(5), cdf.Percentile(95)\n",
    "    print('confidence interval', ci)\n",
    "    VertLine(ci[0])\n",
    "    VertLine(ci[1])\n",
    "\n",
    "    # plot the CDF\n",
    "    thinkplot.Cdf(cdf)\n",
    "    thinkplot.Config(xlabel='estimate',\n",
    "                     ylabel='CDF',\n",
    "                     title='Sampling distribution')\n",
    "\n",
    "    return stderr\n",
    "\n",
    "SimulateSample()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0c9affa4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# My conclusions:\n",
    "\n",
    "# 1) With sample size 10:\n",
    "\n",
    "# standard error 0.762510819389\n",
    "# confidence interval (1.2674054394352277, 3.5377353792673705)\n",
    "\n",
    "# 2) As sample size increases, standard error and the width of\n",
    "#    the CI decrease:\n",
    "\n",
    "# 10      0.90    (1.3, 3.9)\n",
    "# 100     0.21    (1.7, 2.4)\n",
    "# 1000    0.06    (1.9, 2.1)\n",
    "\n",
    "# All three confidence intervals contain the actual value, 2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8d4dc5c",
   "metadata": {},
   "source": [
    "**Exercise:** In games like hockey and soccer, the time between goals is roughly exponential. So you could estimate a team’s goal-scoring rate by observing the number of goals they score in a game. This estimation process is a little different from sampling the time between goals, so let’s see how it works.\n",
    "\n",
    "Write a function that takes a goal-scoring rate, `lam`, in goals per game, and simulates a game by generating the time between goals until the total time exceeds 1 game, then returns the number of goals scored.\n",
    "\n",
    "Write another function that simulates many games, stores the estimates of `lam`, then computes their mean error and RMSE.\n",
    "\n",
    "Is this way of making an estimate biased?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "9dcfda77",
   "metadata": {},
   "outputs": [],
   "source": [
    "def SimulateGame(lam):\n",
    "    \"\"\"Simulates a game and returns the estimated goal-scoring rate.\n",
    "\n",
    "    lam: actual goal scoring rate in goals per game\n",
    "    \"\"\"\n",
    "    goals = 0\n",
    "    t = 0\n",
    "    while True:\n",
    "        time_between_goals = random.expovariate(lam)\n",
    "        t += time_between_goals\n",
    "        if t > 1:\n",
    "            break\n",
    "        goals += 1\n",
    "\n",
    "    # estimated goal-scoring rate is the actual number of goals scored\n",
    "    L = goals\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "2a7a3a99",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 4\n",
      "rmse L 1.4167646240642797\n",
      "mean error L 0.000938\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGwCAYAAABIC3rIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAArSElEQVR4nO3dfVRVZaLH8d8B5YCg+EKCGAqmpRhqSXB96Xpbngmdlult8oVrSdQ4q5Sbxo2MroJraaJmRqbJ2GTqGt9m7qR3mmkoY8TbC4pJVL7kWFmgBKipCCQY7PtHy9OcBBMVtvh8P2vtlWef5+z97J3Rd+2zD8dhWZYlAAAAg3nZPQEAAAC7EUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMF4buydwLaqvr1dJSYnat28vh8Nh93QAAMAlsCxLZ86cUWhoqLy8mnbNhyBqQElJicLCwuyeBgAAuAzFxcW68cYbm/QagqgB7du3l/TDCe3QoYPNswEAAJeioqJCYWFh7v+PNwVB1IDzb5N16NCBIAIAoJW5nNtduKkaAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDx2tg9AbS8pPkbW3yfy2fH27r/n84BAIB/xhUiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGO+aCKIVK1YoPDxcvr6+io2NVX5+fqNjX3nlFd15553q1KmTOnXqJJfLdcH4hx56SA6Hw2MZNWpUcx8GAABopWwPos2bNys5OVnp6ekqKCjQwIEDFRcXp/Ly8gbH5+bmKj4+Xtu3b1deXp7CwsJ099136+jRox7jRo0apW+++ca9bNy4sSUOBwAAtEK2B9HSpUs1depUJSYmKjIyUllZWWrXrp1Wr17d4Pj169dr2rRpGjRokPr27avf/e53qq+vV05Ojsc4p9OpkJAQ99KpU6dG51BTU6OKigqPBQAAmMPWIKqtrdWePXvkcrnc67y8vORyuZSXl3dJ26iurta5c+fUuXNnj/W5ubnq2rWrbrnlFj322GM6ceJEo9vIyMhQYGCgewkLC7u8AwIAAK2SrUF0/Phx1dXVKTg42GN9cHCwSktLL2kbs2bNUmhoqEdUjRo1SuvWrVNOTo4WLVqkHTt2aPTo0aqrq2twG6mpqTp9+rR7KS4uvvyDAgAArU4buydwJRYuXKhNmzYpNzdXvr6+7vWTJk1y/zkqKkoDBgzQTTfdpNzcXI0cOfKC7TidTjmdzhaZMwAAuPbYeoUoKChI3t7eKisr81hfVlamkJCQi752yZIlWrhwod5++20NGDDgomN79eqloKAgff7551c8ZwAAcP2xNYh8fHw0ePBgjxuiz98gPWTIkEZft3jxYs2bN0/Z2dmKjo7+2f0cOXJEJ06cULdu3a7KvAEAwPXF9k+ZJScn65VXXtHatWt14MABPfbYY6qqqlJiYqIkacqUKUpNTXWPX7RokebMmaPVq1crPDxcpaWlKi0tVWVlpSSpsrJSKSkp2rlzp7766ivl5ORo7Nix6t27t+Li4mw5RgAAcG2z/R6iiRMn6tixY0pLS1NpaakGDRqk7Oxs943WRUVF8vL6sdtWrlyp2tpa3X///R7bSU9P19y5c+Xt7a1PPvlEa9eu1alTpxQaGqq7775b8+bN4z4hAADQINuDSJKSkpKUlJTU4HO5ubkej7/66quLbsvPz09vvfXWVZoZAAAwge1vmQEAANiNIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxiOIAACA8QgiAABgPIIIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCAAAGI8gAgAAxrsmgmjFihUKDw+Xr6+vYmNjlZ+f3+jYV155RXfeeac6deqkTp06yeVyXTDesiylpaWpW7du8vPzk8vl0qFDh5r7MAAAQCtlexBt3rxZycnJSk9PV0FBgQYOHKi4uDiVl5c3OD43N1fx8fHavn278vLyFBYWprvvvltHjx51j1m8eLGWLVumrKws7dq1S/7+/oqLi9PZs2db6rAAAEArYnsQLV26VFOnTlViYqIiIyOVlZWldu3aafXq1Q2OX79+vaZNm6ZBgwapb9+++t3vfqf6+nrl5ORI+uHqUGZmpmbPnq2xY8dqwIABWrdunUpKSrR169YGt1lTU6OKigqPBQAAmMPWIKqtrdWePXvkcrnc67y8vORyuZSXl3dJ26iurta5c+fUuXNnSdLhw4dVWlrqsc3AwEDFxsY2us2MjAwFBga6l7CwsCs4KgAA0NrYGkTHjx9XXV2dgoODPdYHBwertLT0krYxa9YshYaGugPo/Ouass3U1FSdPn3avRQXFzf1UAAAQCvWxu4JXImFCxdq06ZNys3Nla+v72Vvx+l0yul0XsWZAQCA1sTWK0RBQUHy9vZWWVmZx/qysjKFhIRc9LVLlizRwoUL9fbbb2vAgAHu9edfdznbBAAAZrI1iHx8fDR48GD3DdGS3DdIDxkypNHXLV68WPPmzVN2draio6M9nouIiFBISIjHNisqKrRr166LbhMAAJjL9rfMkpOTlZCQoOjoaMXExCgzM1NVVVVKTEyUJE2ZMkXdu3dXRkaGJGnRokVKS0vThg0bFB4e7r4vKCAgQAEBAXI4HJo5c6bmz5+vPn36KCIiQnPmzFFoaKjGjRtn12ECAIBrmO1BNHHiRB07dkxpaWkqLS3VoEGDlJ2d7b4puqioSF5eP17IWrlypWpra3X//fd7bCc9PV1z586VJD311FOqqqrSb37zG506dUrDhw9Xdnb2Fd1nBAAArl8Oy7IsuydxramoqFBgYKBOnz6tDh062D2dqy5p/sYW3+fy2fG27v+ncwAAXH+u5P/ftv9iRgAAALsRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA47WxewKAHZLmb2zxfS6fHd/i+wQAXBquEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxnexCtWLFC4eHh8vX1VWxsrPLz8xsdu2/fPv3qV79SeHi4HA6HMjMzLxgzd+5cORwOj6Vv377NeAQAAKC1szWINm/erOTkZKWnp6ugoEADBw5UXFycysvLGxxfXV2tXr16aeHChQoJCWl0u/3799c333zjXt57773mOgQAAHAdsDWIli5dqqlTpyoxMVGRkZHKyspSu3bttHr16gbH33HHHXruuec0adIkOZ3ORrfbpk0bhYSEuJegoKCLzqOmpkYVFRUeCwAAMIdtQVRbW6s9e/bI5XL9OBkvL7lcLuXl5V3Rtg8dOqTQ0FD16tVLkydPVlFR0UXHZ2RkKDAw0L2EhYVd0f4BAEDr0qQgSktLU3V1tfvxyZMnL3vHx48fV11dnYKDgz3WBwcHq7S09LK3GxsbqzVr1ig7O1srV67U4cOHdeedd+rMmTONviY1NVWnT592L8XFxZe9fwAA0Po0KYieffZZVVZWuh/37NlTX3755VWf1JUYPXq0xo8frwEDBiguLk5vvvmmTp06pT/84Q+NvsbpdKpDhw4eCwAAMEeTgsiyrIs+boqgoCB5e3urrKzMY31ZWdlFb5huqo4dO+rmm2/W559/ftW2CQAAri+23UPk4+OjwYMHKycnx72uvr5eOTk5GjJkyFXbT2Vlpb744gt169btqm0TAABcX9o0ZbDD4dCZM2fk6+sry7LkcDhUWVl5waeyLvUtp+TkZCUkJCg6OloxMTHKzMxUVVWVEhMTJUlTpkxR9+7dlZGRIemHG7H379/v/vPRo0dVWFiogIAA9e7dW5L05JNPasyYMerZs6dKSkqUnp4ub29vxcfHN+VQAQCAQZoURJZl6eabb/Z4fNttt3k8djgcqquru6TtTZw4UceOHVNaWppKS0s1aNAgZWdnu2+0LioqkpfXjxexSkpKPPa3ZMkSLVmyRCNGjFBubq4k6ciRI4qPj9eJEyd0ww03aPjw4dq5c6duuOGGphwqAAAwSJOCaPv27Vd9AklJSUpKSmrwufORc154ePjP3re0adOmqzU1AABgiCYF0YgRI5prHgAAALax/bvMAAAA7NakK0Te3t6XNO5S7yECAAC4FjT5puqePXsqISHB4+ZmAACA1qxJQZSfn69XX31VL774oiIiIvTwww9r8uTJ6tSpU3PNDwAAoNk16R6i6OhorVy5Ut98842Sk5O1ZcsW3XjjjZo0aZK2bdvWXHMEAABoVpd1U7Wvr68eeOAB5eTkaO/evSovL9eoUaP07bffXu35AQAANLsmvWX2z44cOaI1a9ZozZo1qq6uVkpKCl+KCgAAWqUmBVFtba22bNmiV199Ve+++65Gjx6tzMxMjR49+pI/gQYAAHCtaVIQdevWTe3bt1dCQoJefvllde3aVZJUVVXlMY4rRQAAoDVpUhCdPHlSJ0+e1Lx58zR//vwLnm/qd5kBAABcC2z/LjMTJc3f2OL7XD47vsX3CQBAa9GkIBo+fLiWLFmiP//5z6qtrdXIkSOVnp4uPz+/5pofAABAs2vSx+4XLFigZ555RgEBAerevbtefPFFTZ8+vbnmBgAA0CKaFETr1q3Tyy+/rLfeektbt27VG2+8ofXr16u+vr655gcAANDsmhRERUVF+uUvf+l+7HK55HA4VFJSctUnBgAA0FKaFETff/+9fH19Pda1bdtW586du6qTAgAAaElN/rb7hx56SE6n073u7NmzevTRR+Xv7+9e9/rrr1+9GQIAADSzJgVRQkLCBeseeOCBqzYZAAAAOzQpiF577bXmmgcAAIBtLuvb7gEAAK4nBBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwnu1BtGLFCoWHh8vX11exsbHKz89vdOy+ffv0q1/9SuHh4XI4HMrMzLzibQIAANgaRJs3b1ZycrLS09NVUFCggQMHKi4uTuXl5Q2Or66uVq9evbRw4UKFhIRclW0CAADYGkRLly7V1KlTlZiYqMjISGVlZaldu3ZavXp1g+PvuOMOPffcc5o0aZKcTudV2SYAAIBtQVRbW6s9e/bI5XL9OBkvL7lcLuXl5bXoNmtqalRRUeGxAAAAc9gWRMePH1ddXZ2Cg4M91gcHB6u0tLRFt5mRkaHAwED3EhYWdln7BwAArZPtN1VfC1JTU3X69Gn3UlxcbPeUAABAC2pj146DgoLk7e2tsrIyj/VlZWWN3jDdXNt0Op2N3pMEAACuf7ZdIfLx8dHgwYOVk5PjXldfX6+cnBwNGTLkmtkmAAC4/tl2hUiSkpOTlZCQoOjoaMXExCgzM1NVVVVKTEyUJE2ZMkXdu3dXRkaGpB9umt6/f7/7z0ePHlVhYaECAgLUu3fvS9omAADAT9kaRBMnTtSxY8eUlpam0tJSDRo0SNnZ2e6boouKiuTl9eNFrJKSEt12223ux0uWLNGSJUs0YsQI5ebmXtI2AQAAfsrWIJKkpKQkJSUlNfjc+cg5Lzw8XJZlXdE2AQAAfopPmQEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA49n+KTPAREnzN9qy3+Wz423ZLwBc67hCBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMN41EUQrVqxQeHi4fH19FRsbq/z8/IuO/+Mf/6i+ffvK19dXUVFRevPNNz2ef+ihh+RwODyWUaNGNechAACAVsz2INq8ebOSk5OVnp6ugoICDRw4UHFxcSovL29w/AcffKD4+Hg98sgj+uijjzRu3DiNGzdOe/fu9Rg3atQoffPNN+5l48aNLXE4AACgFbI9iJYuXaqpU6cqMTFRkZGRysrKUrt27bR69eoGx7/44osaNWqUUlJS1K9fP82bN0+33367li9f7jHO6XQqJCTEvXTq1KklDgcAALRCtgZRbW2t9uzZI5fL5V7n5eUll8ulvLy8Bl+Tl5fnMV6S4uLiLhifm5urrl276pZbbtFjjz2mEydONDqPmpoaVVRUeCwAAMActgbR8ePHVVdXp+DgYI/1wcHBKi0tbfA1paWlPzt+1KhRWrdunXJycrRo0SLt2LFDo0ePVl1dXYPbzMjIUGBgoHsJCwu7wiMDAACtSRu7J9AcJk2a5P5zVFSUBgwYoJtuukm5ubkaOXLkBeNTU1OVnJzsflxRUUEUAQBgEFuvEAUFBcnb21tlZWUe68vKyhQSEtLga0JCQpo0XpJ69eqloKAgff755w0+73Q61aFDB48FAACYw9Yg8vHx0eDBg5WTk+NeV19fr5ycHA0ZMqTB1wwZMsRjvCRt27at0fGSdOTIEZ04cULdunW7OhMHAADXFds/ZZacnKxXXnlFa9eu1YEDB/TYY4+pqqpKiYmJkqQpU6YoNTXVPX7GjBnKzs7W888/r88++0xz587Vhx9+qKSkJElSZWWlUlJStHPnTn311VfKycnR2LFj1bt3b8XFxdlyjAAA4Npm+z1EEydO1LFjx5SWlqbS0lINGjRI2dnZ7huni4qK5OX1Y7cNHTpUGzZs0OzZs/XMM8+oT58+2rp1q2699VZJkre3tz755BOtXbtWp06dUmhoqO6++27NmzdPTqfTlmMEAADXNtuDSJKSkpLcV3h+Kjc394J148eP1/jx4xsc7+fnp7feeutqTg8AAFznbH/LDAAAwG4EEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIx3TXyXGYCWlzR/Y4vvc/ns+BbfJwBcCq4QAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADjEUQAAMB4beyeAAAzJc3f2OL7XD47vsX3CaB14AoRAAAwHkEEAACMRxABAADjEUQAAMB4BBEAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjMdXdwAwFl8fAuC8a+IK0YoVKxQeHi5fX1/FxsYqPz//ouP/+Mc/qm/fvvL19VVUVJTefPNNj+cty1JaWpq6desmPz8/uVwuHTp0qDkPAQAAtGK2B9HmzZuVnJys9PR0FRQUaODAgYqLi1N5eXmD4z/44APFx8frkUce0UcffaRx48Zp3Lhx2rt3r3vM4sWLtWzZMmVlZWnXrl3y9/dXXFyczp4921KHBQAAWhHbg2jp0qWaOnWqEhMTFRkZqaysLLVr106rV69ucPyLL76oUaNGKSUlRf369dO8efN0++23a/ny5ZJ+uDqUmZmp2bNna+zYsRowYIDWrVunkpISbd26tQWPDAAAtBa23kNUW1urPXv2KDU11b3Oy8tLLpdLeXl5Db4mLy9PycnJHuvi4uLcsXP48GGVlpbK5XK5nw8MDFRsbKzy8vI0adKkC7ZZU1Ojmpoa9+PTp09LkioqKi772C6m9mx1s2z3Yv75WEzc/7UwB7v3fy3M4Vra/7UwhycX/7HF97/kqfEej6+FOQBXy/n/vizLavqLLRsdPXrUkmR98MEHHutTUlKsmJiYBl/Ttm1ba8OGDR7rVqxYYXXt2tWyLMt6//33LUlWSUmJx5jx48dbEyZMaHCb6enpliQWFhYWFhaW62ApLi5ucpPwKTNJqampHled6uvr9e2336pLly5yOBw2zuxHFRUVCgsLU3FxsTp06GD3dGzBOeAcSJwDiXMgcQ5MP36p4XNgWZbOnDmj0NDQJm/P1iAKCgqSt7e3ysrKPNaXlZUpJCSkwdeEhIRcdPz5f5aVlalbt24eYwYNGtTgNp1Op5xOp8e6jh07NuVQWkyHDh2M/ct/HueAcyBxDiTOgcQ5MP34pQvPQWBg4GVtx9abqn18fDR48GDl5OS419XX1ysnJ0dDhgxp8DVDhgzxGC9J27Ztc4+PiIhQSEiIx5iKigrt2rWr0W0CAACz2f6WWXJyshISEhQdHa2YmBhlZmaqqqpKiYmJkqQpU6aoe/fuysjIkCTNmDFDI0aM0PPPP6977rlHmzZt0ocffqhVq1ZJkhwOh2bOnKn58+erT58+ioiI0Jw5cxQaGqpx48bZdZgAAOAaZnsQTZw4UceOHVNaWppKS0s1aNAgZWdnKzg4WJJUVFQkL68fL2QNHTpUGzZs0OzZs/XMM8+oT58+2rp1q2699Vb3mKeeekpVVVX6zW9+o1OnTmn48OHKzs6Wr69vix/f1eJ0OpWenn7BW3sm4RxwDiTOgcQ5kDgHph+/dPXPgcOyLuezaQAAANcP238xIwAAgN0IIgAAYDyCCAAAGI8gAgAAxiOIWokVK1YoPDxcvr6+io2NVX5+vt1TajEZGRm644471L59e3Xt2lXjxo3TwYMH7Z6WbRYuXOj+9RImOXr0qB544AF16dJFfn5+ioqK0ocffmj3tFpMXV2d5syZo4iICPn5+emmm27SvHnzLu87m1qJ//u//9OYMWMUGhoqh8NxwRd0W5altLQ0devWTX5+fnK5XDp06JA9k20mFzsH586d06xZsxQVFSV/f3+FhoZqypQpKikpsW/CzeDn/h78s0cffVQOh0OZmZlN3g9B1Aps3rxZycnJSk9PV0FBgQYOHKi4uDiVl5fbPbUWsWPHDk2fPl07d+7Utm3bdO7cOd19992qqqqye2otbvfu3frtb3+rAQMG2D2VFnXy5EkNGzZMbdu21d/+9jft379fzz//vDp16mT31FrMokWLtHLlSi1fvlwHDhzQokWLtHjxYr300kt2T63ZVFVVaeDAgVqxYkWDzy9evFjLli1TVlaWdu3aJX9/f8XFxens2bMtPNPmc7FzUF1drYKCAs2ZM0cFBQV6/fXXdfDgQd177702zLT5/Nzfg/O2bNminTt3XtbXdkiSrV/uiksTExNjTZ8+3f24rq7OCg0NtTIyMmyclX3Ky8stSdaOHTvsnkqLOnPmjNWnTx9r27Zt1ogRI6wZM2bYPaUWM2vWLGv48OF2T8NW99xzj/Xwww97rLvvvvusyZMn2zSjliXJ2rJli/txfX29FRISYj333HPudadOnbKcTqe1ceNGG2bY/H56DhqSn59vSbK+/vrrlplUC2vsHBw5csTq3r27tXfvXqtnz57WCy+80ORtc4XoGldbW6s9e/bI5XK513l5ecnlcikvL8/Gmdnn9OnTkqTOnTvbPJOWNX36dN1zzz0efxdM8ec//1nR0dEaP368unbtqttuu02vvPKK3dNqUUOHDlVOTo7+8Y9/SJI+/vhjvffeexo9erTNM7PH4cOHVVpa6vHfQ2BgoGJjY4392Sj98PPR4XBcs9/H2Rzq6+v14IMPKiUlRf3797/s7dj+m6pxccePH1ddXZ37N3efFxwcrM8++8ymWdmnvr5eM2fO1LBhwzx+O/n1btOmTSooKNDu3bvtnootvvzyS61cuVLJycl65plntHv3bj3++OPy8fFRQkKC3dNrEU8//bQqKirUt29feXt7q66uTs8++6wmT55s99RsUVpaKkkN/mw8/5xpzp49q1mzZik+Pt6oL3xdtGiR2rRpo8cff/yKtkMQoVWZPn269u7dq/fee8/uqbSY4uJizZgxQ9u2bWvVXz9zJerr6xUdHa0FCxZIkm677Tbt3btXWVlZxgTRH/7wB61fv14bNmxQ//79VVhYqJkzZyo0NNSYc4DGnTt3ThMmTJBlWVq5cqXd02kxe/bs0YsvvqiCggI5HI4r2hZvmV3jgoKC5O3trbKyMo/1ZWVlCgkJsWlW9khKStJf/vIXbd++XTfeeKPd02kxe/bsUXl5uW6//Xa1adNGbdq00Y4dO7Rs2TK1adNGdXV1dk+x2XXr1k2RkZEe6/r166eioiKbZtTyUlJS9PTTT2vSpEmKiorSgw8+qCeeeML9xdemOf/zj5+NP8bQ119/rW3bthl1dejdd99VeXm5evTo4f75+PXXX+u//uu/FB4e3qRtEUTXOB8fHw0ePFg5OTnudfX19crJydGQIUNsnFnLsSxLSUlJ2rJli/7+978rIiLC7im1qJEjR+rTTz9VYWGhe4mOjtbkyZNVWFgob29vu6fY7IYNG3bBr1r4xz/+oZ49e9o0o5ZXXV3t8UXXkuTt7a36+nqbZmSviIgIhYSEePxsrKio0K5du4z52Sj9GEOHDh3SO++8oy5dutg9pRb14IMP6pNPPvH4+RgaGqqUlBS99dZbTdoWb5m1AsnJyUpISFB0dLRiYmKUmZmpqqoqJSYm2j21FjF9+nRt2LBB//u//6v27du77w8IDAyUn5+fzbNrfu3bt7/gfil/f3916dLFmPuonnjiCQ0dOlQLFizQhAkTlJ+fr1WrVmnVqlV2T63FjBkzRs8++6x69Oih/v3766OPPtLSpUv18MMP2z21ZlNZWanPP//c/fjw4cMqLCxU586d1aNHD82cOVPz589Xnz59FBERoTlz5ig0NFTjxo2zb9JX2cXOQbdu3XT//feroKBAf/nLX1RXV+f++di5c2f5+PjYNe2r6uf+Hvw0Atu2bauQkBDdcsstTdvRlX4EDi3jpZdesnr06GH5+PhYMTEx1s6dO+2eUouR1ODy2muv2T0125j2sXvLsqw33njDuvXWWy2n02n17dvXWrVqld1TalEVFRXWjBkzrB49eli+vr5Wr169rP/+7/+2ampq7J5as9m+fXuD/+0nJCRYlvXDR+/nzJljBQcHW06n0xo5cqR18OBBeyd9lV3sHBw+fLjRn4/bt2+3e+pXzc/9Pfipy/3YvcOyruNfcwoAAHAJuIcIAAAYjyACAADGI4gAAIDxCCIAAGA8gggAABiPIAIAAMYjiAAAgPEIIgAAYDyCCECrEx4erszMTLun0axyc3PlcDh06tQpu6cCGIEgAnBVlJaWasaMGerdu7d8fX0VHBysYcOGaeXKlaqurrZ7egBwUXy5K4Ar9uWXX2rYsGHq2LGjFixYoKioKDmdTn366adatWqVunfvrnvvvdfuabY4y7JUV1enNm34UQtc67hCBOCKTZs2TW3atNGHH36oCRMmqF+/furVq5fGjh2rv/71rxozZox7bFFRkcaOHauAgAB16NBBEyZMUFlZmfv5L774QmPHjlVwcLACAgJ0xx136J133ml035Zlae7cuerRo4ecTqdCQ0P1+OOPNzr+448/1l133aX27durQ4cOGjx4sD788EP38++//77+7d/+Te3atVOnTp0UFxenkydPSpJqamr0+OOPq2vXrvL19dXw4cO1e/du92vPv831t7/9TYMHD5bT6dR7772n+vp6ZWRkKCIiQn5+fho4cKD+53/+x2Neb775pm6++Wb5+fnprrvu0ldffXXJ5x/AlSOIAFyREydO6O2339b06dPl7+/f4BiHwyFJqq+v19ixY/Xtt99qx44d2rZtm7788ktNnDjRPbayslK//OUvlZOTo48++kijRo3SmDFjVFRU1OC2//SnP+mFF17Qb3/7Wx06dEhbt25VVFRUo/OdPHmybrzxRu3evVt79uzR008/rbZt20qSCgsLNXLkSEVGRiovL0/vvfeexowZo7q6OknSU089pT/96U9au3atCgoK1Lt3b8XFxenbb7/12MfTTz+thQsX6sCBAxowYIAyMjK0bt06ZWVlad++fXriiSf0wAMPaMeOHZKk4uJi3XfffRozZowKCwv161//Wk8//fQl/hsAcFVYAHAFdu7caUmyXn/9dY/1Xbp0sfz9/S1/f3/rqaeesizLst5++23L29vbKioqco/bt2+fJcnKz89vdB/9+/e3XnrpJffjnj17Wi+88IJlWZb1/PPPWzfffLNVW1t7SfNt3769tWbNmgafi4+Pt4YNG9bgc5WVlVbbtm2t9evXu9fV1tZaoaGh1uLFiy3Lsqzt27dbkqytW7e6x5w9e9Zq166d9cEHH3hs75FHHrHi4+Mty7Ks1NRUKzIy0uP5WbNmWZKskydPXtJxAbgyXCEC0Czy8/NVWFio/v37q6amRpJ04MABhYWFKSwszD0uMjJSHTt21IEDByT9cIXoySefVL9+/dSxY0cFBATowIEDjV4hGj9+vL777jv16tVLU6dO1ZYtW/T99983Oq/k5GT9+te/lsvl0sKFC/XFF1+4nzt/haghX3zxhc6dO6dhw4a517Vt21YxMTHuuZ8XHR3t/vPnn3+u6upq/eIXv1BAQIB7WbdunXvfBw4cUGxsrMc2hgwZ0ugxALj6CCIAV6R3795yOBw6ePCgx/pevXqpd+/e8vPza9L2nnzySW3ZskULFizQu+++q8LCQkVFRam2trbB8WFhYTp48KBefvll+fn5adq0afrXf/1XnTt3rsHxc+fO1b59+3TPPffo73//uyIjI7VlyxZJavJcG/PPbx1WVlZKkv7617+qsLDQvezfv/+C+4gA2IcgAnBFunTpol/84hdavny5qqqqLjq2X79+Ki4uVnFxsXvd/v37derUKUVGRkr64abmhx56SP/+7/+uqKgohYSE/OwNxn5+fhozZoyWLVum3Nxc5eXl6dNPP210/M0336wnnnhCb7/9tu677z699tprkqQBAwYoJyenwdfcdNNN8vHx0fvvv+9ed+7cOe3evds994ZERkbK6XSqqKhIvXv39ljOXynr16+f8vPzPV63c+fOix4zgKuLIAJwxV5++WV9//33io6O1ubNm3XgwAEdPHhQv//97/XZZ5/J29tbkuRyuRQVFaXJkyeroKBA+fn5mjJlikaMGOF+m6lPnz56/fXXVVhYqI8//lj/8R//ofr6+kb3vWbNGr366qvau3evvvzyS/3+97+Xn5+fevbsecHY7777TklJScrNzdXXX3+t999/X7t371a/fv0kSampqdq9e7emTZumTz75RJ999plWrlyp48ePy9/fX4899phSUlKUnZ2t/fv3a+rUqaqurtYjjzzS6Pzat2+vJ598Uk888YTWrl2rL774QgUFBXrppZe0du1aSdKjjz6qQ4cOKSUlRQcPHtSGDRu0Zs2ay/3XAeBy2H0TE4DrQ0lJiZWUlGRFRERYbdu2tQICAqyYmBjrueees6qqqtzjvv76a+vee++1/P39rfbt21vjx4+3SktL3c8fPnzYuuuuuyw/Pz8rLCzMWr58uTVixAhrxowZ7jH/fFP1li1brNjYWKtDhw6Wv7+/9S//8i/WO++80+Aca2pqrEmTJllhYWGWj4+PFRoaaiUlJVnfffede0xubq41dOhQy+l0Wh07drTi4uLcNzZ/99131n/+539aQUFBltPptIYNG+ZxM/j5m6p/eiN0fX29lZmZad1yyy1W27ZtrRtuuMGKi4uzduzY4R7zxhtvWL1797acTqd15513WqtXr+amaqAFOSzLsuyOMgAAADvxlhkAADAeQQQAAIxHEAEAAOMRRAAAwHgEEQAAMB5BBAAAjEcQAQAA4xFEAADAeAQRAAAwHkEEAACMRxABAADj/T/OG114igRF6gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# The following function simulates many games, then uses the\n",
    "# number of goals scored as an estimate of the true long-term\n",
    "# goal-scoring rate.\n",
    "\n",
    "def Estimate6(lam=2, m=1000000):\n",
    "\n",
    "    estimates = []\n",
    "    for i in range(m):\n",
    "        L = SimulateGame(lam)\n",
    "        estimates.append(L)\n",
    "\n",
    "    print('Experiment 4')\n",
    "    print('rmse L', RMSE(estimates, lam))\n",
    "    print('mean error L', MeanError(estimates, lam))\n",
    "    \n",
    "    pmf = thinkstats2.Pmf(estimates)\n",
    "    thinkplot.Hist(pmf)\n",
    "    thinkplot.Config(xlabel='Goals scored', ylabel='PMF')\n",
    "    \n",
    "Estimate6()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "1ec0bff3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# My conclusions:\n",
    "\n",
    "# 1) RMSE for this way of estimating lambda is 1.4\n",
    "\n",
    "# 2) The mean error is small and decreases with m, so this estimator\n",
    "#    appears to be unbiased.\n",
    "\n",
    "# One note: If the time between goals is exponential, the distribution\n",
    "# of goals scored in a game is Poisson.\n",
    "\n",
    "# See https://en.wikipedia.org/wiki/Poisson_distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f4fbe65",
   "metadata": {},
   "source": [
    "**Exercise:**  In this chapter we used $\\bar{x}$ and median to estimate µ, and found that $\\bar{x}$ yields lower MSE. Also, we used $S^2$ and $S_{n-1}^2$ to estimate σ, and found that $S^2$ is biased and $S_{n-1}^2$ unbiased.\n",
    "Run similar experiments to see if $\\bar{x}$ and median are biased estimates of µ. Also check whether $S^2$ or $S_{n-1}^2$ yields a lower MSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "52ac2dbc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 1\n",
      "mean error xbar -0.0004525331787192879\n",
      "mean error median -0.0006871767159844678\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "def Estimate4(n=7, iters=100000):\n",
    "    \"\"\"Mean error for xbar and median as estimators of population mean.\n",
    "\n",
    "    n: sample size\n",
    "    iters: number of iterations\n",
    "    \"\"\"\n",
    "    mu = 0\n",
    "    sigma = 1\n",
    "\n",
    "    means = []\n",
    "    medians = []\n",
    "    for _ in range(iters):\n",
    "        xs = [random.gauss(mu, sigma) for i in range(n)]\n",
    "        xbar = np.mean(xs)\n",
    "        median = np.median(xs)\n",
    "        means.append(xbar)\n",
    "        medians.append(median)\n",
    "\n",
    "    print('Experiment 1')\n",
    "    print('mean error xbar', MeanError(means, mu))\n",
    "    print('mean error median', MeanError(medians, mu))\n",
    "    \n",
    "Estimate4()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "bbffceac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Experiment 2\n",
      "RMSE biased 0.5128622081175854\n",
      "RMSE unbiased 0.5744124496387355\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "def Estimate5(n=7, iters=100000):\n",
    "    \"\"\"RMSE for biased and unbiased estimators of population variance.\n",
    "\n",
    "    n: sample size\n",
    "    iters: number of iterations\n",
    "    \"\"\"\n",
    "    mu = 0\n",
    "    sigma = 1\n",
    "\n",
    "    estimates1 = []\n",
    "    estimates2 = []\n",
    "    for _ in range(iters):\n",
    "        xs = [random.gauss(mu, sigma) for i in range(n)]\n",
    "        biased = np.var(xs)\n",
    "        unbiased = np.var(xs, ddof=1)\n",
    "        estimates1.append(biased)\n",
    "        estimates2.append(unbiased)\n",
    "\n",
    "    print('Experiment 2')\n",
    "    print('RMSE biased', RMSE(estimates1, sigma**2))\n",
    "    print('RMSE unbiased', RMSE(estimates2, sigma**2))\n",
    "\n",
    "Estimate5()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "bac4114e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# My conclusions:\n",
    "\n",
    "# 1) xbar and median yield lower mean error as m increases, so neither\n",
    "# one is obviously biased, as far as we can tell from the experiment.\n",
    "\n",
    "# 2) The biased estimator of variance yields lower RMSE than the unbiased\n",
    "# estimator, by about 10%.  And the difference holds up as m increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37084486",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
